Actual source code: ispai.c


  2: /*
  3:    3/99 Modified by Stephen Barnard to support SPAI version 3.0
  4: */

  6: /*
  7:       Provides an interface to the SPAI Sparse Approximate Inverse Preconditioner
  8:    Code written by Stephen Barnard.

 10:       Note: there is some BAD memory bleeding below!

 12:       This code needs work

 14:    1) get rid of all memory bleeding
 15:    2) fix PETSc/interface so that it gets if the matrix is symmetric from the matrix
 16:       rather than having the sp flag for PC_SPAI
 17:    3) fix to set the block size based on the matrix block size

 19: */
 20: #define PETSC_SKIP_COMPLEX /* since spai uses I which conflicts with some complex implementations */

 22: #include <petsc/private/pcimpl.h>
 23: #include <../src/ksp/pc/impls/spai/petscspai.h>

 25: /*
 26:     These are the SPAI include files
 27: */
 28: EXTERN_C_BEGIN
 29: #define SPAI_USE_MPI /* required for setting SPAI_Comm correctly in basics.h */
 30: #include <spai.h>
 31: #include <matrix.h>
 32: EXTERN_C_END

 34: extern PetscErrorCode ConvertMatToMatrix(MPI_Comm,Mat,Mat,matrix**);
 35: extern PetscErrorCode ConvertMatrixToMat(MPI_Comm,matrix*,Mat*);
 36: extern PetscErrorCode ConvertVectorToVec(MPI_Comm,vector *v,Vec *Pv);
 37: extern PetscErrorCode MM_to_PETSC(char*,char*,char*);

 39: typedef struct {

 41:   matrix *B;                /* matrix in SPAI format */
 42:   matrix *BT;               /* transpose of matrix in SPAI format */
 43:   matrix *M;                /* the approximate inverse in SPAI format */

 45:   Mat PM;                   /* the approximate inverse PETSc format */

 47:   double epsilon;           /* tolerance */
 48:   int    nbsteps;           /* max number of "improvement" steps per line */
 49:   int    max;               /* max dimensions of is_I, q, etc. */
 50:   int    maxnew;            /* max number of new entries per step */
 51:   int    block_size;        /* constant block size */
 52:   int    cache_size;        /* one of (1,2,3,4,5,6) indicting size of cache */
 53:   int    verbose;           /* SPAI prints timing and statistics */

 55:   int      sp;              /* symmetric nonzero pattern */
 56:   MPI_Comm comm_spai;     /* communicator to be used with spai */
 57: } PC_SPAI;

 59: /**********************************************************************/

 61: static PetscErrorCode PCSetUp_SPAI(PC pc)
 62: {
 63:   PC_SPAI *ispai = (PC_SPAI*)pc->data;
 64:   Mat      AT;

 66:   init_SPAI();

 68:   if (ispai->sp) {
 69:     ConvertMatToMatrix(ispai->comm_spai,pc->pmat,pc->pmat,&ispai->B);
 70:   } else {
 71:     /* Use the transpose to get the column nonzero structure. */
 72:     MatTranspose(pc->pmat,MAT_INITIAL_MATRIX,&AT);
 73:     ConvertMatToMatrix(ispai->comm_spai,pc->pmat,AT,&ispai->B);
 74:     MatDestroy(&AT);
 75:   }

 77:   /* Destroy the transpose */
 78:   /* Don't know how to do it. PETSc developers? */

 80:   /* construct SPAI preconditioner */
 81:   /* FILE *messages */     /* file for warning messages */
 82:   /* double epsilon */     /* tolerance */
 83:   /* int nbsteps */        /* max number of "improvement" steps per line */
 84:   /* int max */            /* max dimensions of is_I, q, etc. */
 85:   /* int maxnew */         /* max number of new entries per step */
 86:   /* int block_size */     /* block_size == 1 specifies scalar elments
 87:                               block_size == n specifies nxn constant-block elements
 88:                               block_size == 0 specifies variable-block elements */
 89:   /* int cache_size */     /* one of (1,2,3,4,5,6) indicting size of cache. cache_size == 0 indicates no caching */
 90:   /* int    verbose    */  /* verbose == 0 specifies that SPAI is silent
 91:                               verbose == 1 prints timing and matrix statistics */

 93:   PetscCall(bspai(ispai->B,&ispai->M,
 94:                 stdout,
 95:                 ispai->epsilon,
 96:                 ispai->nbsteps,
 97:                 ispai->max,
 98:                 ispai->maxnew,
 99:                 ispai->block_size,
100:                 ispai->cache_size,
101:                 ispai->verbose));

103:   ConvertMatrixToMat(PetscObjectComm((PetscObject)pc),ispai->M,&ispai->PM);

105:   /* free the SPAI matrices */
106:   sp_free_matrix(ispai->B);
107:   sp_free_matrix(ispai->M);
108:   return 0;
109: }

111: /**********************************************************************/

113: static PetscErrorCode PCApply_SPAI(PC pc,Vec xx,Vec y)
114: {
115:   PC_SPAI *ispai = (PC_SPAI*)pc->data;

117:   /* Now using PETSc's multiply */
118:   MatMult(ispai->PM,xx,y);
119:   return 0;
120: }

122: static PetscErrorCode PCMatApply_SPAI(PC pc,Mat X,Mat Y)
123: {
124:   PC_SPAI *ispai = (PC_SPAI*)pc->data;

126:   /* Now using PETSc's multiply */
127:   MatMatMult(ispai->PM,X,MAT_REUSE_MATRIX,PETSC_DEFAULT,&Y);
128:   return 0;
129: }

131: /**********************************************************************/

133: static PetscErrorCode PCDestroy_SPAI(PC pc)
134: {
135:   PC_SPAI *ispai = (PC_SPAI*)pc->data;

137:   MatDestroy(&ispai->PM);
138:   MPI_Comm_free(&(ispai->comm_spai));
139:   PetscFree(pc->data);
140:   return 0;
141: }

143: /**********************************************************************/

145: static PetscErrorCode PCView_SPAI(PC pc,PetscViewer viewer)
146: {
147:   PC_SPAI   *ispai = (PC_SPAI*)pc->data;
148:   PetscBool  iascii;

150:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
151:   if (iascii) {
152:     PetscViewerASCIIPrintf(viewer,"    epsilon %g\n",   (double)ispai->epsilon);
153:     PetscViewerASCIIPrintf(viewer,"    nbsteps %d\n",   ispai->nbsteps);
154:     PetscViewerASCIIPrintf(viewer,"    max %d\n",       ispai->max);
155:     PetscViewerASCIIPrintf(viewer,"    maxnew %d\n",    ispai->maxnew);
156:     PetscViewerASCIIPrintf(viewer,"    block_size %d\n",ispai->block_size);
157:     PetscViewerASCIIPrintf(viewer,"    cache_size %d\n",ispai->cache_size);
158:     PetscViewerASCIIPrintf(viewer,"    verbose %d\n",   ispai->verbose);
159:     PetscViewerASCIIPrintf(viewer,"    sp %d\n",        ispai->sp);
160:   }
161:   return 0;
162: }

164: static PetscErrorCode  PCSPAISetEpsilon_SPAI(PC pc,double epsilon1)
165: {
166:   PC_SPAI *ispai = (PC_SPAI*)pc->data;

168:   ispai->epsilon = epsilon1;
169:   return 0;
170: }

172: /**********************************************************************/

174: static PetscErrorCode  PCSPAISetNBSteps_SPAI(PC pc,int nbsteps1)
175: {
176:   PC_SPAI *ispai = (PC_SPAI*)pc->data;

178:   ispai->nbsteps = nbsteps1;
179:   return 0;
180: }

182: /**********************************************************************/

184: /* added 1/7/99 g.h. */
185: static PetscErrorCode  PCSPAISetMax_SPAI(PC pc,int max1)
186: {
187:   PC_SPAI *ispai = (PC_SPAI*)pc->data;

189:   ispai->max = max1;
190:   return 0;
191: }

193: /**********************************************************************/

195: static PetscErrorCode  PCSPAISetMaxNew_SPAI(PC pc,int maxnew1)
196: {
197:   PC_SPAI *ispai = (PC_SPAI*)pc->data;

199:   ispai->maxnew = maxnew1;
200:   return 0;
201: }

203: /**********************************************************************/

205: static PetscErrorCode  PCSPAISetBlockSize_SPAI(PC pc,int block_size1)
206: {
207:   PC_SPAI *ispai = (PC_SPAI*)pc->data;

209:   ispai->block_size = block_size1;
210:   return 0;
211: }

213: /**********************************************************************/

215: static PetscErrorCode  PCSPAISetCacheSize_SPAI(PC pc,int cache_size)
216: {
217:   PC_SPAI *ispai = (PC_SPAI*)pc->data;

219:   ispai->cache_size = cache_size;
220:   return 0;
221: }

223: /**********************************************************************/

225: static PetscErrorCode  PCSPAISetVerbose_SPAI(PC pc,int verbose)
226: {
227:   PC_SPAI *ispai = (PC_SPAI*)pc->data;

229:   ispai->verbose = verbose;
230:   return 0;
231: }

233: /**********************************************************************/

235: static PetscErrorCode  PCSPAISetSp_SPAI(PC pc,int sp)
236: {
237:   PC_SPAI *ispai = (PC_SPAI*)pc->data;

239:   ispai->sp = sp;
240:   return 0;
241: }

243: /* -------------------------------------------------------------------*/

245: /*@
246:   PCSPAISetEpsilon -- Set the tolerance for the SPAI preconditioner

248:   Input Parameters:
249: + pc - the preconditioner
250: - eps - epsilon (default .4)

252:   Notes:
253:     Espilon must be between 0 and 1. It controls the
254:                  quality of the approximation of M to the inverse of
255:                  A. Higher values of epsilon lead to more work, more
256:                  fill, and usually better preconditioners. In many
257:                  cases the best choice of epsilon is the one that
258:                  divides the total solution time equally between the
259:                  preconditioner and the solver.

261:   Level: intermediate

263: .seealso: PCSPAI, PCSetType()
264:   @*/
265: PetscErrorCode  PCSPAISetEpsilon(PC pc,double epsilon1)
266: {
267:   PetscTryMethod(pc,"PCSPAISetEpsilon_C",(PC,double),(pc,epsilon1));
268:   return 0;
269: }

271: /**********************************************************************/

273: /*@
274:   PCSPAISetNBSteps - set maximum number of improvement steps per row in
275:         the SPAI preconditioner

277:   Input Parameters:
278: + pc - the preconditioner
279: - n - number of steps (default 5)

281:   Notes:
282:     SPAI constructs to approximation to every column of
283:                  the exact inverse of A in a series of improvement
284:                  steps. The quality of the approximation is determined
285:                  by epsilon. If an approximation achieving an accuracy
286:                  of epsilon is not obtained after ns steps, SPAI simply
287:                  uses the best approximation constructed so far.

289:   Level: intermediate

291: .seealso: PCSPAI, PCSetType(), PCSPAISetMaxNew()
292: @*/
293: PetscErrorCode  PCSPAISetNBSteps(PC pc,int nbsteps1)
294: {
295:   PetscTryMethod(pc,"PCSPAISetNBSteps_C",(PC,int),(pc,nbsteps1));
296:   return 0;
297: }

299: /**********************************************************************/

301: /* added 1/7/99 g.h. */
302: /*@
303:   PCSPAISetMax - set the size of various working buffers in
304:         the SPAI preconditioner

306:   Input Parameters:
307: + pc - the preconditioner
308: - n - size (default is 5000)

310:   Level: intermediate

312: .seealso: PCSPAI, PCSetType()
313: @*/
314: PetscErrorCode  PCSPAISetMax(PC pc,int max1)
315: {
316:   PetscTryMethod(pc,"PCSPAISetMax_C",(PC,int),(pc,max1));
317:   return 0;
318: }

320: /**********************************************************************/

322: /*@
323:   PCSPAISetMaxNew - set maximum number of new nonzero candidates per step
324:    in SPAI preconditioner

326:   Input Parameters:
327: + pc - the preconditioner
328: - n - maximum number (default 5)

330:   Level: intermediate

332: .seealso: PCSPAI, PCSetType(), PCSPAISetNBSteps()
333: @*/
334: PetscErrorCode  PCSPAISetMaxNew(PC pc,int maxnew1)
335: {
336:   PetscTryMethod(pc,"PCSPAISetMaxNew_C",(PC,int),(pc,maxnew1));
337:   return 0;
338: }

340: /**********************************************************************/

342: /*@
343:   PCSPAISetBlockSize - set the block size for the SPAI preconditioner

345:   Input Parameters:
346: + pc - the preconditioner
347: - n - block size (default 1)

349:   Notes:
350:     A block
351:                  size of 1 treats A as a matrix of scalar elements. A
352:                  block size of s > 1 treats A as a matrix of sxs
353:                  blocks. A block size of 0 treats A as a matrix with
354:                  variable sized blocks, which are determined by
355:                  searching for dense square diagonal blocks in A.
356:                  This can be very effective for finite-element
357:                  matrices.

359:                  SPAI will convert A to block form, use a block
360:                  version of the preconditioner algorithm, and then
361:                  convert the result back to scalar form.

363:                  In many cases the a block-size parameter other than 1
364:                  can lead to very significant improvement in
365:                  performance.

367:   Level: intermediate

369: .seealso: PCSPAI, PCSetType()
370: @*/
371: PetscErrorCode  PCSPAISetBlockSize(PC pc,int block_size1)
372: {
373:   PetscTryMethod(pc,"PCSPAISetBlockSize_C",(PC,int),(pc,block_size1));
374:   return 0;
375: }

377: /**********************************************************************/

379: /*@
380:   PCSPAISetCacheSize - specify cache size in the SPAI preconditioner

382:   Input Parameters:
383: + pc - the preconditioner
384: - n -  cache size {0,1,2,3,4,5} (default 5)

386:   Notes:
387:     SPAI uses a hash table to cache messages and avoid
388:                  redundant communication. If suggest always using
389:                  5. This parameter is irrelevant in the serial
390:                  version.

392:   Level: intermediate

394: .seealso: PCSPAI, PCSetType()
395: @*/
396: PetscErrorCode  PCSPAISetCacheSize(PC pc,int cache_size)
397: {
398:   PetscTryMethod(pc,"PCSPAISetCacheSize_C",(PC,int),(pc,cache_size));
399:   return 0;
400: }

402: /**********************************************************************/

404: /*@
405:   PCSPAISetVerbose - verbosity level for the SPAI preconditioner

407:   Input Parameters:
408: + pc - the preconditioner
409: - n - level (default 1)

411:   Notes:
412:     print parameters, timings and matrix statistics

414:   Level: intermediate

416: .seealso: PCSPAI, PCSetType()
417: @*/
418: PetscErrorCode  PCSPAISetVerbose(PC pc,int verbose)
419: {
420:   PetscTryMethod(pc,"PCSPAISetVerbose_C",(PC,int),(pc,verbose));
421:   return 0;
422: }

424: /**********************************************************************/

426: /*@
427:   PCSPAISetSp - specify a symmetric matrix sparsity pattern in the SPAI preconditioner

429:   Input Parameters:
430: + pc - the preconditioner
431: - n - 0 or 1

433:   Notes:
434:     If A has a symmetric nonzero pattern use -sp 1 to
435:                  improve performance by eliminating some communication
436:                  in the parallel version. Even if A does not have a
437:                  symmetric nonzero pattern -sp 1 may well lead to good
438:                  results, but the code will not follow the published
439:                  SPAI algorithm exactly.

441:   Level: intermediate

443: .seealso: PCSPAI, PCSetType()
444: @*/
445: PetscErrorCode  PCSPAISetSp(PC pc,int sp)
446: {
447:   PetscTryMethod(pc,"PCSPAISetSp_C",(PC,int),(pc,sp));
448:   return 0;
449: }

451: /**********************************************************************/

453: /**********************************************************************/

455: static PetscErrorCode PCSetFromOptions_SPAI(PetscOptionItems *PetscOptionsObject,PC pc)
456: {
457:   PC_SPAI        *ispai = (PC_SPAI*)pc->data;
458:   int            nbsteps1,max1,maxnew1,block_size1,cache_size,verbose,sp;
459:   double         epsilon1;
460:   PetscBool      flg;

462:   PetscOptionsHead(PetscOptionsObject,"SPAI options");
463:   PetscOptionsReal("-pc_spai_epsilon","","PCSPAISetEpsilon",ispai->epsilon,&epsilon1,&flg);
464:   if (flg) {
465:     PCSPAISetEpsilon(pc,epsilon1);
466:   }
467:   PetscOptionsInt("-pc_spai_nbsteps","","PCSPAISetNBSteps",ispai->nbsteps,&nbsteps1,&flg);
468:   if (flg) {
469:     PCSPAISetNBSteps(pc,nbsteps1);
470:   }
471:   /* added 1/7/99 g.h. */
472:   PetscOptionsInt("-pc_spai_max","","PCSPAISetMax",ispai->max,&max1,&flg);
473:   if (flg) {
474:     PCSPAISetMax(pc,max1);
475:   }
476:   PetscOptionsInt("-pc_spai_maxnew","","PCSPAISetMaxNew",ispai->maxnew,&maxnew1,&flg);
477:   if (flg) {
478:     PCSPAISetMaxNew(pc,maxnew1);
479:   }
480:   PetscOptionsInt("-pc_spai_block_size","","PCSPAISetBlockSize",ispai->block_size,&block_size1,&flg);
481:   if (flg) {
482:     PCSPAISetBlockSize(pc,block_size1);
483:   }
484:   PetscOptionsInt("-pc_spai_cache_size","","PCSPAISetCacheSize",ispai->cache_size,&cache_size,&flg);
485:   if (flg) {
486:     PCSPAISetCacheSize(pc,cache_size);
487:   }
488:   PetscOptionsInt("-pc_spai_verbose","","PCSPAISetVerbose",ispai->verbose,&verbose,&flg);
489:   if (flg) {
490:     PCSPAISetVerbose(pc,verbose);
491:   }
492:   PetscOptionsInt("-pc_spai_sp","","PCSPAISetSp",ispai->sp,&sp,&flg);
493:   if (flg) {
494:     PCSPAISetSp(pc,sp);
495:   }
496:   PetscOptionsTail();
497:   return 0;
498: }

500: /**********************************************************************/

502: /*MC
503:    PCSPAI - Use the Sparse Approximate Inverse method of Grote and Barnard
504:      as a preconditioner (SIAM J. Sci. Comput.; vol 18, nr 3)

506:    Options Database Keys:
507: +  -pc_spai_epsilon <eps> - set tolerance
508: .  -pc_spai_nbstep <n> - set nbsteps
509: .  -pc_spai_max <m> - set max
510: .  -pc_spai_max_new <m> - set maxnew
511: .  -pc_spai_block_size <n> - set block size
512: .  -pc_spai_cache_size <n> - set cache size
513: .  -pc_spai_sp <m> - set sp
514: -  -pc_spai_set_verbose <true,false> - verbose output

516:    Notes:
517:     This only works with AIJ matrices.

519:    Level: beginner

521: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC,
522:     PCSPAISetEpsilon(), PCSPAISetMax(), PCSPAISetMaxNew(), PCSPAISetBlockSize(),
523:     PCSPAISetVerbose(), PCSPAISetSp()
524: M*/

526: PETSC_EXTERN PetscErrorCode PCCreate_SPAI(PC pc)
527: {
528:   PC_SPAI *ispai;

530:   PetscNewLog(pc,&ispai);
531:   pc->data = ispai;

533:   pc->ops->destroy         = PCDestroy_SPAI;
534:   pc->ops->apply           = PCApply_SPAI;
535:   pc->ops->matapply        = PCMatApply_SPAI;
536:   pc->ops->applyrichardson = 0;
537:   pc->ops->setup           = PCSetUp_SPAI;
538:   pc->ops->view            = PCView_SPAI;
539:   pc->ops->setfromoptions  = PCSetFromOptions_SPAI;

541:   ispai->epsilon    = .4;
542:   ispai->nbsteps    = 5;
543:   ispai->max        = 5000;
544:   ispai->maxnew     = 5;
545:   ispai->block_size = 1;
546:   ispai->cache_size = 5;
547:   ispai->verbose    = 0;

549:   ispai->sp = 1;
550:   MPI_Comm_dup(PetscObjectComm((PetscObject)pc),&(ispai->comm_spai));

552:   PetscObjectComposeFunction((PetscObject)pc,"PCSPAISetEpsilon_C",PCSPAISetEpsilon_SPAI);
553:   PetscObjectComposeFunction((PetscObject)pc,"PCSPAISetNBSteps_C",PCSPAISetNBSteps_SPAI);
554:   PetscObjectComposeFunction((PetscObject)pc,"PCSPAISetMax_C",PCSPAISetMax_SPAI);
555:   PetscObjectComposeFunction((PetscObject)pc,"PCSPAISetMaxNew_C",PCSPAISetMaxNew_SPAI);
556:   PetscObjectComposeFunction((PetscObject)pc,"PCSPAISetBlockSize_C",PCSPAISetBlockSize_SPAI);
557:   PetscObjectComposeFunction((PetscObject)pc,"PCSPAISetCacheSize_C",PCSPAISetCacheSize_SPAI);
558:   PetscObjectComposeFunction((PetscObject)pc,"PCSPAISetVerbose_C",PCSPAISetVerbose_SPAI);
559:   PetscObjectComposeFunction((PetscObject)pc,"PCSPAISetSp_C",PCSPAISetSp_SPAI);
560:   return 0;
561: }

563: /**********************************************************************/

565: /*
566:    Converts from a PETSc matrix to an SPAI matrix
567: */
568: PetscErrorCode ConvertMatToMatrix(MPI_Comm comm, Mat A,Mat AT,matrix **B)
569: {
570:   matrix                  *M;
571:   int                     i,j,col;
572:   int                     row_indx;
573:   int                     len,pe,local_indx,start_indx;
574:   int                     *mapping;
575:   const int               *cols;
576:   const double            *vals;
577:   int                     n,mnl,nnl,nz,rstart,rend;
578:   PetscMPIInt             size,rank;
579:   struct compressed_lines *rows;

581:   MPI_Comm_size(comm,&size);
582:   MPI_Comm_rank(comm,&rank);
583:   MatGetSize(A,&n,&n);
584:   MatGetLocalSize(A,&mnl,&nnl);

586:   /*
587:     not sure why a barrier is required. commenting out
588:   MPI_Barrier(comm);
589:   */

591:   M = new_matrix((SPAI_Comm)comm);

593:   M->n              = n;
594:   M->bs             = 1;
595:   M->max_block_size = 1;

597:   M->mnls          = (int*)malloc(sizeof(int)*size);
598:   M->start_indices = (int*)malloc(sizeof(int)*size);
599:   M->pe            = (int*)malloc(sizeof(int)*n);
600:   M->block_sizes   = (int*)malloc(sizeof(int)*n);
601:   for (i=0; i<n; i++) M->block_sizes[i] = 1;

603:   MPI_Allgather(&mnl,1,MPI_INT,M->mnls,1,MPI_INT,comm);

605:   M->start_indices[0] = 0;
606:   for (i=1; i<size; i++) M->start_indices[i] = M->start_indices[i-1] + M->mnls[i-1];

608:   M->mnl            = M->mnls[M->myid];
609:   M->my_start_index = M->start_indices[M->myid];

611:   for (i=0; i<size; i++) {
612:     start_indx = M->start_indices[i];
613:     for (j=0; j<M->mnls[i]; j++) M->pe[start_indx+j] = i;
614:   }

616:   if (AT) {
617:     M->lines = new_compressed_lines(M->mnls[rank],1);
618:   } else {
619:     M->lines = new_compressed_lines(M->mnls[rank],0);
620:   }

622:   rows = M->lines;

624:   /* Determine the mapping from global indices to pointers */
625:   PetscMalloc1(M->n,&mapping);
626:   pe         = 0;
627:   local_indx = 0;
628:   for (i=0; i<M->n; i++) {
629:     if (local_indx >= M->mnls[pe]) {
630:       pe++;
631:       local_indx = 0;
632:     }
633:     mapping[i] = local_indx + M->start_indices[pe];
634:     local_indx++;
635:   }

637:   /*********************************************************/
638:   /************** Set up the row structure *****************/
639:   /*********************************************************/

641:   MatGetOwnershipRange(A,&rstart,&rend);
642:   for (i=rstart; i<rend; i++) {
643:     row_indx = i - rstart;
644:     MatGetRow(A,i,&nz,&cols,&vals);
645:     /* allocate buffers */
646:     rows->ptrs[row_indx] = (int*)malloc(nz*sizeof(int));
647:     rows->A[row_indx]    = (double*)malloc(nz*sizeof(double));
648:     /* copy the matrix */
649:     for (j=0; j<nz; j++) {
650:       col = cols[j];
651:       len = rows->len[row_indx]++;

653:       rows->ptrs[row_indx][len] = mapping[col];
654:       rows->A[row_indx][len]    = vals[j];
655:     }
656:     rows->slen[row_indx] = rows->len[row_indx];

658:     MatRestoreRow(A,i,&nz,&cols,&vals);
659:   }

661:   /************************************************************/
662:   /************** Set up the column structure *****************/
663:   /*********************************************************/

665:   if (AT) {

667:     for (i=rstart; i<rend; i++) {
668:       row_indx = i - rstart;
669:       MatGetRow(AT,i,&nz,&cols,&vals);
670:       /* allocate buffers */
671:       rows->rptrs[row_indx] = (int*)malloc(nz*sizeof(int));
672:       /* copy the matrix (i.e., the structure) */
673:       for (j=0; j<nz; j++) {
674:         col = cols[j];
675:         len = rows->rlen[row_indx]++;

677:         rows->rptrs[row_indx][len] = mapping[col];
678:       }
679:       MatRestoreRow(AT,i,&nz,&cols,&vals);
680:     }
681:   }

683:   PetscFree(mapping);

685:   order_pointers(M);
686:   M->maxnz = calc_maxnz(M);
687:   *B       = M;
688:   return 0;
689: }

691: /**********************************************************************/

693: /*
694:    Converts from an SPAI matrix B  to a PETSc matrix PB.
695:    This assumes that the SPAI matrix B is stored in
696:    COMPRESSED-ROW format.
697: */
698: PetscErrorCode ConvertMatrixToMat(MPI_Comm comm,matrix *B,Mat *PB)
699: {
700:   PetscMPIInt    size,rank;
701:   int            m,n,M,N;
702:   int            d_nz,o_nz;
703:   int            *d_nnz,*o_nnz;
704:   int            i,k,global_row,global_col,first_diag_col,last_diag_col;
705:   PetscScalar    val;

707:   MPI_Comm_size(comm,&size);
708:   MPI_Comm_rank(comm,&rank);

710:   m    = n = B->mnls[rank];
711:   d_nz = o_nz = 0;

713:   /* Determine preallocation for MatCreateAIJ */
714:   PetscMalloc1(m,&d_nnz);
715:   PetscMalloc1(m,&o_nnz);
716:   for (i=0; i<m; i++) d_nnz[i] = o_nnz[i] = 0;
717:   first_diag_col = B->start_indices[rank];
718:   last_diag_col  = first_diag_col + B->mnls[rank];
719:   for (i=0; i<B->mnls[rank]; i++) {
720:     for (k=0; k<B->lines->len[i]; k++) {
721:       global_col = B->lines->ptrs[i][k];
722:       if ((global_col >= first_diag_col) && (global_col < last_diag_col)) d_nnz[i]++;
723:       else o_nnz[i]++;
724:     }
725:   }

727:   M = N = B->n;
728:   /* Here we only know how to create AIJ format */
729:   MatCreate(comm,PB);
730:   MatSetSizes(*PB,m,n,M,N);
731:   MatSetType(*PB,MATAIJ);
732:   MatSeqAIJSetPreallocation(*PB,d_nz,d_nnz);
733:   MatMPIAIJSetPreallocation(*PB,d_nz,d_nnz,o_nz,o_nnz);

735:   for (i=0; i<B->mnls[rank]; i++) {
736:     global_row = B->start_indices[rank]+i;
737:     for (k=0; k<B->lines->len[i]; k++) {
738:       global_col = B->lines->ptrs[i][k];

740:       val  = B->lines->A[i][k];
741:       MatSetValues(*PB,1,&global_row,1,&global_col,&val,ADD_VALUES);
742:     }
743:   }

745:   PetscFree(d_nnz);
746:   PetscFree(o_nnz);

748:   MatAssemblyBegin(*PB,MAT_FINAL_ASSEMBLY);
749:   MatAssemblyEnd(*PB,MAT_FINAL_ASSEMBLY);
750:   return 0;
751: }

753: /**********************************************************************/

755: /*
756:    Converts from an SPAI vector v  to a PETSc vec Pv.
757: */
758: PetscErrorCode ConvertVectorToVec(MPI_Comm comm,vector *v,Vec *Pv)
759: {
760:   PetscMPIInt size,rank;
761:   int         m,M,i,*mnls,*start_indices,*global_indices;

763:   MPI_Comm_size(comm,&size);
764:   MPI_Comm_rank(comm,&rank);

766:   m = v->mnl;
767:   M = v->n;

769:   VecCreateMPI(comm,m,M,Pv);

771:   PetscMalloc1(size,&mnls);
772:   MPI_Allgather(&v->mnl,1,MPI_INT,mnls,1,MPI_INT,comm);

774:   PetscMalloc1(size,&start_indices);

776:   start_indices[0] = 0;
777:   for (i=1; i<size; i++) start_indices[i] = start_indices[i-1] +mnls[i-1];

779:   PetscMalloc1(v->mnl,&global_indices);
780:   for (i=0; i<v->mnl; i++) global_indices[i] = start_indices[rank] + i;

782:   PetscFree(mnls);
783:   PetscFree(start_indices);

785:   VecSetValues(*Pv,v->mnl,global_indices,v->v,INSERT_VALUES);
786:   VecAssemblyBegin(*Pv);
787:   VecAssemblyEnd(*Pv);

789:   PetscFree(global_indices);
790:   return 0;
791: }