Actual source code: ispai.c

petsc-3.11.4 2019-09-28
Report Typos and Errors

  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;
 65:   Mat            AT;

 68:   init_SPAI();

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

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

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

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

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

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

113: /**********************************************************************/

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

121:   /* Now using PETSc's multiply */
122:   MatMult(ispai->PM,xx,y);
123:   return(0);
124: }

126: /**********************************************************************/

128: static PetscErrorCode PCDestroy_SPAI(PC pc)
129: {
131:   PC_SPAI        *ispai = (PC_SPAI*)pc->data;

134:   MatDestroy(&ispai->PM);
135:   MPI_Comm_free(&(ispai->comm_spai));
136:   PetscFree(pc->data);
137:   return(0);
138: }

140: /**********************************************************************/

142: static PetscErrorCode PCView_SPAI(PC pc,PetscViewer viewer)
143: {
144:   PC_SPAI        *ispai = (PC_SPAI*)pc->data;
146:   PetscBool      iascii;

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

163: static PetscErrorCode  PCSPAISetEpsilon_SPAI(PC pc,double epsilon1)
164: {
165:   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;

179:   ispai->nbsteps = nbsteps1;
180:   return(0);
181: }

183: /**********************************************************************/

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

191:   ispai->max = max1;
192:   return(0);
193: }

195: /**********************************************************************/

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

202:   ispai->maxnew = maxnew1;
203:   return(0);
204: }

206: /**********************************************************************/

208: static PetscErrorCode  PCSPAISetBlockSize_SPAI(PC pc,int block_size1)
209: {
210:   PC_SPAI *ispai = (PC_SPAI*)pc->data;

213:   ispai->block_size = block_size1;
214:   return(0);
215: }

217: /**********************************************************************/

219: static PetscErrorCode  PCSPAISetCacheSize_SPAI(PC pc,int cache_size)
220: {
221:   PC_SPAI *ispai = (PC_SPAI*)pc->data;

224:   ispai->cache_size = cache_size;
225:   return(0);
226: }

228: /**********************************************************************/

230: static PetscErrorCode  PCSPAISetVerbose_SPAI(PC pc,int verbose)
231: {
232:   PC_SPAI *ispai = (PC_SPAI*)pc->data;

235:   ispai->verbose = verbose;
236:   return(0);
237: }

239: /**********************************************************************/

241: static PetscErrorCode  PCSPAISetSp_SPAI(PC pc,int sp)
242: {
243:   PC_SPAI *ispai = (PC_SPAI*)pc->data;

246:   ispai->sp = sp;
247:   return(0);
248: }

250: /* -------------------------------------------------------------------*/

252: /*@
253:   PCSPAISetEpsilon -- Set the tolerance for the SPAI preconditioner

255:   Input Parameters:
256: + pc - the preconditioner
257: - eps - epsilon (default .4)

259:   Notes:
260:     Espilon must be between 0 and 1. It controls the
261:                  quality of the approximation of M to the inverse of
262:                  A. Higher values of epsilon lead to more work, more
263:                  fill, and usually better preconditioners. In many
264:                  cases the best choice of epsilon is the one that
265:                  divides the total solution time equally between the
266:                  preconditioner and the solver.

268:   Level: intermediate

270: .seealso: PCSPAI, PCSetType()
271:   @*/
272: PetscErrorCode  PCSPAISetEpsilon(PC pc,double epsilon1)
273: {

277:   PetscTryMethod(pc,"PCSPAISetEpsilon_C",(PC,double),(pc,epsilon1));
278:   return(0);
279: }

281: /**********************************************************************/

283: /*@
284:   PCSPAISetNBSteps - set maximum number of improvement steps per row in
285:         the SPAI preconditioner

287:   Input Parameters:
288: + pc - the preconditioner
289: - n - number of steps (default 5)

291:   Notes:
292:     SPAI constructs to approximation to every column of
293:                  the exact inverse of A in a series of improvement
294:                  steps. The quality of the approximation is determined
295:                  by epsilon. If an approximation achieving an accuracy
296:                  of epsilon is not obtained after ns steps, SPAI simply
297:                  uses the best approximation constructed so far.

299:   Level: intermediate

301: .seealso: PCSPAI, PCSetType(), PCSPAISetMaxNew()
302: @*/
303: PetscErrorCode  PCSPAISetNBSteps(PC pc,int nbsteps1)
304: {

308:   PetscTryMethod(pc,"PCSPAISetNBSteps_C",(PC,int),(pc,nbsteps1));
309:   return(0);
310: }

312: /**********************************************************************/

314: /* added 1/7/99 g.h. */
315: /*@
316:   PCSPAISetMax - set the size of various working buffers in
317:         the SPAI preconditioner

319:   Input Parameters:
320: + pc - the preconditioner
321: - n - size (default is 5000)

323:   Level: intermediate

325: .seealso: PCSPAI, PCSetType()
326: @*/
327: PetscErrorCode  PCSPAISetMax(PC pc,int max1)
328: {

332:   PetscTryMethod(pc,"PCSPAISetMax_C",(PC,int),(pc,max1));
333:   return(0);
334: }

336: /**********************************************************************/

338: /*@
339:   PCSPAISetMaxNew - set maximum number of new nonzero candidates per step
340:    in SPAI preconditioner

342:   Input Parameters:
343: + pc - the preconditioner
344: - n - maximum number (default 5)

346:   Level: intermediate

348: .seealso: PCSPAI, PCSetType(), PCSPAISetNBSteps()
349: @*/
350: PetscErrorCode  PCSPAISetMaxNew(PC pc,int maxnew1)
351: {

355:   PetscTryMethod(pc,"PCSPAISetMaxNew_C",(PC,int),(pc,maxnew1));
356:   return(0);
357: }

359: /**********************************************************************/

361: /*@
362:   PCSPAISetBlockSize - set the block size for the SPAI preconditioner

364:   Input Parameters:
365: + pc - the preconditioner
366: - n - block size (default 1)

368:   Notes:
369:     A block
370:                  size of 1 treats A as a matrix of scalar elements. A
371:                  block size of s > 1 treats A as a matrix of sxs
372:                  blocks. A block size of 0 treats A as a matrix with
373:                  variable sized blocks, which are determined by
374:                  searching for dense square diagonal blocks in A.
375:                  This can be very effective for finite-element
376:                  matrices.

378:                  SPAI will convert A to block form, use a block
379:                  version of the preconditioner algorithm, and then
380:                  convert the result back to scalar form.

382:                  In many cases the a block-size parameter other than 1
383:                  can lead to very significant improvement in
384:                  performance.


387:   Level: intermediate

389: .seealso: PCSPAI, PCSetType()
390: @*/
391: PetscErrorCode  PCSPAISetBlockSize(PC pc,int block_size1)
392: {

396:   PetscTryMethod(pc,"PCSPAISetBlockSize_C",(PC,int),(pc,block_size1));
397:   return(0);
398: }

400: /**********************************************************************/

402: /*@
403:   PCSPAISetCacheSize - specify cache size in the SPAI preconditioner

405:   Input Parameters:
406: + pc - the preconditioner
407: - n -  cache size {0,1,2,3,4,5} (default 5)

409:   Notes:
410:     SPAI uses a hash table to cache messages and avoid
411:                  redundant communication. If suggest always using
412:                  5. This parameter is irrelevant in the serial
413:                  version.

415:   Level: intermediate

417: .seealso: PCSPAI, PCSetType()
418: @*/
419: PetscErrorCode  PCSPAISetCacheSize(PC pc,int cache_size)
420: {

424:   PetscTryMethod(pc,"PCSPAISetCacheSize_C",(PC,int),(pc,cache_size));
425:   return(0);
426: }

428: /**********************************************************************/

430: /*@
431:   PCSPAISetVerbose - verbosity level for the SPAI preconditioner

433:   Input Parameters:
434: + pc - the preconditioner
435: - n - level (default 1)

437:   Notes:
438:     print parameters, timings and matrix statistics

440:   Level: intermediate

442: .seealso: PCSPAI, PCSetType()
443: @*/
444: PetscErrorCode  PCSPAISetVerbose(PC pc,int verbose)
445: {

449:   PetscTryMethod(pc,"PCSPAISetVerbose_C",(PC,int),(pc,verbose));
450:   return(0);
451: }

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

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

458:   Input Parameters:
459: + pc - the preconditioner
460: - n - 0 or 1

462:   Notes:
463:     If A has a symmetric nonzero pattern use -sp 1 to
464:                  improve performance by eliminating some communication
465:                  in the parallel version. Even if A does not have a
466:                  symmetric nonzero pattern -sp 1 may well lead to good
467:                  results, but the code will not follow the published
468:                  SPAI algorithm exactly.


471:   Level: intermediate

473: .seealso: PCSPAI, PCSetType()
474: @*/
475: PetscErrorCode  PCSPAISetSp(PC pc,int sp)
476: {

480:   PetscTryMethod(pc,"PCSPAISetSp_C",(PC,int),(pc,sp));
481:   return(0);
482: }

484: /**********************************************************************/

486: /**********************************************************************/

488: static PetscErrorCode PCSetFromOptions_SPAI(PetscOptionItems *PetscOptionsObject,PC pc)
489: {
490:   PC_SPAI        *ispai = (PC_SPAI*)pc->data;
492:   int            nbsteps1,max1,maxnew1,block_size1,cache_size,verbose,sp;
493:   double         epsilon1;
494:   PetscBool      flg;

497:   PetscOptionsHead(PetscOptionsObject,"SPAI options");
498:   PetscOptionsReal("-pc_spai_epsilon","","PCSPAISetEpsilon",ispai->epsilon,&epsilon1,&flg);
499:   if (flg) {
500:     PCSPAISetEpsilon(pc,epsilon1);
501:   }
502:   PetscOptionsInt("-pc_spai_nbsteps","","PCSPAISetNBSteps",ispai->nbsteps,&nbsteps1,&flg);
503:   if (flg) {
504:     PCSPAISetNBSteps(pc,nbsteps1);
505:   }
506:   /* added 1/7/99 g.h. */
507:   PetscOptionsInt("-pc_spai_max","","PCSPAISetMax",ispai->max,&max1,&flg);
508:   if (flg) {
509:     PCSPAISetMax(pc,max1);
510:   }
511:   PetscOptionsInt("-pc_spai_maxnew","","PCSPAISetMaxNew",ispai->maxnew,&maxnew1,&flg);
512:   if (flg) {
513:     PCSPAISetMaxNew(pc,maxnew1);
514:   }
515:   PetscOptionsInt("-pc_spai_block_size","","PCSPAISetBlockSize",ispai->block_size,&block_size1,&flg);
516:   if (flg) {
517:     PCSPAISetBlockSize(pc,block_size1);
518:   }
519:   PetscOptionsInt("-pc_spai_cache_size","","PCSPAISetCacheSize",ispai->cache_size,&cache_size,&flg);
520:   if (flg) {
521:     PCSPAISetCacheSize(pc,cache_size);
522:   }
523:   PetscOptionsInt("-pc_spai_verbose","","PCSPAISetVerbose",ispai->verbose,&verbose,&flg);
524:   if (flg) {
525:     PCSPAISetVerbose(pc,verbose);
526:   }
527:   PetscOptionsInt("-pc_spai_sp","","PCSPAISetSp",ispai->sp,&sp,&flg);
528:   if (flg) {
529:     PCSPAISetSp(pc,sp);
530:   }
531:   PetscOptionsTail();
532:   return(0);
533: }

535: /**********************************************************************/

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

541:    Options Database Keys:
542: +  -pc_spai_epsilon <eps> - set tolerance
543: .  -pc_spai_nbstep <n> - set nbsteps
544: .  -pc_spai_max <m> - set max
545: .  -pc_spai_max_new <m> - set maxnew
546: .  -pc_spai_block_size <n> - set block size
547: .  -pc_spai_cache_size <n> - set cache size
548: .  -pc_spai_sp <m> - set sp
549: -  -pc_spai_set_verbose <true,false> - verbose output

551:    Notes:
552:     This only works with AIJ matrices.

554:    Level: beginner

556:    Concepts: approximate inverse

558: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC,
559:     PCSPAISetEpsilon(), PCSPAISetMax(), PCSPAISetMaxNew(), PCSPAISetBlockSize(),
560:     PCSPAISetVerbose(), PCSPAISetSp()
561: M*/

563: PETSC_EXTERN PetscErrorCode PCCreate_SPAI(PC pc)
564: {
565:   PC_SPAI        *ispai;

569:   PetscNewLog(pc,&ispai);
570:   pc->data = ispai;

572:   pc->ops->destroy         = PCDestroy_SPAI;
573:   pc->ops->apply           = PCApply_SPAI;
574:   pc->ops->applyrichardson = 0;
575:   pc->ops->setup           = PCSetUp_SPAI;
576:   pc->ops->view            = PCView_SPAI;
577:   pc->ops->setfromoptions  = PCSetFromOptions_SPAI;

579:   ispai->epsilon    = .4;
580:   ispai->nbsteps    = 5;
581:   ispai->max        = 5000;
582:   ispai->maxnew     = 5;
583:   ispai->block_size = 1;
584:   ispai->cache_size = 5;
585:   ispai->verbose    = 0;

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

590:   PetscObjectComposeFunction((PetscObject)pc,"PCSPAISetEpsilon_C",PCSPAISetEpsilon_SPAI);
591:   PetscObjectComposeFunction((PetscObject)pc,"PCSPAISetNBSteps_C",PCSPAISetNBSteps_SPAI);
592:   PetscObjectComposeFunction((PetscObject)pc,"PCSPAISetMax_C",PCSPAISetMax_SPAI);
593:   PetscObjectComposeFunction((PetscObject)pc,"PCSPAISetMaxNew_C",PCSPAISetMaxNew_SPAI);
594:   PetscObjectComposeFunction((PetscObject)pc,"PCSPAISetBlockSize_C",PCSPAISetBlockSize_SPAI);
595:   PetscObjectComposeFunction((PetscObject)pc,"PCSPAISetCacheSize_C",PCSPAISetCacheSize_SPAI);
596:   PetscObjectComposeFunction((PetscObject)pc,"PCSPAISetVerbose_C",PCSPAISetVerbose_SPAI);
597:   PetscObjectComposeFunction((PetscObject)pc,"PCSPAISetSp_C",PCSPAISetSp_SPAI);
598:   return(0);
599: }

601: /**********************************************************************/

603: /*
604:    Converts from a PETSc matrix to an SPAI matrix
605: */
606: PetscErrorCode ConvertMatToMatrix(MPI_Comm comm, Mat A,Mat AT,matrix **B)
607: {
608:   matrix                  *M;
609:   int                     i,j,col;
610:   int                     row_indx;
611:   int                     len,pe,local_indx,start_indx;
612:   int                     *mapping;
613:   PetscErrorCode          ierr;
614:   const int               *cols;
615:   const double            *vals;
616:   int                     n,mnl,nnl,nz,rstart,rend;
617:   PetscMPIInt             size,rank;
618:   struct compressed_lines *rows;

621:   MPI_Comm_size(comm,&size);
622:   MPI_Comm_rank(comm,&rank);
623:   MatGetSize(A,&n,&n);
624:   MatGetLocalSize(A,&mnl,&nnl);

626:   /*
627:     not sure why a barrier is required. commenting out
628:   MPI_Barrier(comm);
629:   */

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

633:   M->n              = n;
634:   M->bs             = 1;
635:   M->max_block_size = 1;

637:   M->mnls          = (int*)malloc(sizeof(int)*size);
638:   M->start_indices = (int*)malloc(sizeof(int)*size);
639:   M->pe            = (int*)malloc(sizeof(int)*n);
640:   M->block_sizes   = (int*)malloc(sizeof(int)*n);
641:   for (i=0; i<n; i++) M->block_sizes[i] = 1;

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

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

648:   M->mnl            = M->mnls[M->myid];
649:   M->my_start_index = M->start_indices[M->myid];

651:   for (i=0; i<size; i++) {
652:     start_indx = M->start_indices[i];
653:     for (j=0; j<M->mnls[i]; j++) M->pe[start_indx+j] = i;
654:   }

656:   if (AT) {
657:     M->lines = new_compressed_lines(M->mnls[rank],1);
658:   } else {
659:     M->lines = new_compressed_lines(M->mnls[rank],0);
660:   }

662:   rows = M->lines;

664:   /* Determine the mapping from global indices to pointers */
665:   PetscMalloc1(M->n,&mapping);
666:   pe         = 0;
667:   local_indx = 0;
668:   for (i=0; i<M->n; i++) {
669:     if (local_indx >= M->mnls[pe]) {
670:       pe++;
671:       local_indx = 0;
672:     }
673:     mapping[i] = local_indx + M->start_indices[pe];
674:     local_indx++;
675:   }

677:   /*********************************************************/
678:   /************** Set up the row structure *****************/
679:   /*********************************************************/

681:   MatGetOwnershipRange(A,&rstart,&rend);
682:   for (i=rstart; i<rend; i++) {
683:     row_indx = i - rstart;
684:     MatGetRow(A,i,&nz,&cols,&vals);
685:     /* allocate buffers */
686:     rows->ptrs[row_indx] = (int*)malloc(nz*sizeof(int));
687:     rows->A[row_indx]    = (double*)malloc(nz*sizeof(double));
688:     /* copy the matrix */
689:     for (j=0; j<nz; j++) {
690:       col = cols[j];
691:       len = rows->len[row_indx]++;

693:       rows->ptrs[row_indx][len] = mapping[col];
694:       rows->A[row_indx][len]    = vals[j];
695:     }
696:     rows->slen[row_indx] = rows->len[row_indx];

698:     MatRestoreRow(A,i,&nz,&cols,&vals);
699:   }


702:   /************************************************************/
703:   /************** Set up the column structure *****************/
704:   /*********************************************************/

706:   if (AT) {

708:     for (i=rstart; i<rend; i++) {
709:       row_indx = i - rstart;
710:       MatGetRow(AT,i,&nz,&cols,&vals);
711:       /* allocate buffers */
712:       rows->rptrs[row_indx] = (int*)malloc(nz*sizeof(int));
713:       /* copy the matrix (i.e., the structure) */
714:       for (j=0; j<nz; j++) {
715:         col = cols[j];
716:         len = rows->rlen[row_indx]++;

718:         rows->rptrs[row_indx][len] = mapping[col];
719:       }
720:       MatRestoreRow(AT,i,&nz,&cols,&vals);
721:     }
722:   }

724:   PetscFree(mapping);

726:   order_pointers(M);
727:   M->maxnz = calc_maxnz(M);
728:   *B       = M;
729:   return(0);
730: }

732: /**********************************************************************/

734: /*
735:    Converts from an SPAI matrix B  to a PETSc matrix PB.
736:    This assumes that the SPAI matrix B is stored in
737:    COMPRESSED-ROW format.
738: */
739: PetscErrorCode ConvertMatrixToMat(MPI_Comm comm,matrix *B,Mat *PB)
740: {
741:   PetscMPIInt    size,rank;
743:   int            m,n,M,N;
744:   int            d_nz,o_nz;
745:   int            *d_nnz,*o_nnz;
746:   int            i,k,global_row,global_col,first_diag_col,last_diag_col;
747:   PetscScalar    val;

750:   MPI_Comm_size(comm,&size);
751:   MPI_Comm_rank(comm,&rank);

753:   m    = n = B->mnls[rank];
754:   d_nz = o_nz = 0;

756:   /* Determine preallocation for MatCreateMPIAIJ */
757:   PetscMalloc1(m,&d_nnz);
758:   PetscMalloc1(m,&o_nnz);
759:   for (i=0; i<m; i++) d_nnz[i] = o_nnz[i] = 0;
760:   first_diag_col = B->start_indices[rank];
761:   last_diag_col  = first_diag_col + B->mnls[rank];
762:   for (i=0; i<B->mnls[rank]; i++) {
763:     for (k=0; k<B->lines->len[i]; k++) {
764:       global_col = B->lines->ptrs[i][k];
765:       if ((global_col >= first_diag_col) && (global_col < last_diag_col)) d_nnz[i]++;
766:       else o_nnz[i]++;
767:     }
768:   }

770:   M = N = B->n;
771:   /* Here we only know how to create AIJ format */
772:   MatCreate(comm,PB);
773:   MatSetSizes(*PB,m,n,M,N);
774:   MatSetType(*PB,MATAIJ);
775:   MatSeqAIJSetPreallocation(*PB,d_nz,d_nnz);
776:   MatMPIAIJSetPreallocation(*PB,d_nz,d_nnz,o_nz,o_nnz);

778:   for (i=0; i<B->mnls[rank]; i++) {
779:     global_row = B->start_indices[rank]+i;
780:     for (k=0; k<B->lines->len[i]; k++) {
781:       global_col = B->lines->ptrs[i][k];

783:       val  = B->lines->A[i][k];
784:       MatSetValues(*PB,1,&global_row,1,&global_col,&val,ADD_VALUES);
785:     }
786:   }

788:   PetscFree(d_nnz);
789:   PetscFree(o_nnz);

791:   MatAssemblyBegin(*PB,MAT_FINAL_ASSEMBLY);
792:   MatAssemblyEnd(*PB,MAT_FINAL_ASSEMBLY);
793:   return(0);
794: }

796: /**********************************************************************/

798: /*
799:    Converts from an SPAI vector v  to a PETSc vec Pv.
800: */
801: PetscErrorCode ConvertVectorToVec(MPI_Comm comm,vector *v,Vec *Pv)
802: {
804:   PetscMPIInt    size,rank;
805:   int            m,M,i,*mnls,*start_indices,*global_indices;

808:   MPI_Comm_size(comm,&size);
809:   MPI_Comm_rank(comm,&rank);

811:   m = v->mnl;
812:   M = v->n;


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

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

820:   PetscMalloc1(size,&start_indices);

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

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

828:   PetscFree(mnls);
829:   PetscFree(start_indices);

831:   VecSetValues(*Pv,v->mnl,global_indices,v->v,INSERT_VALUES);
832:   VecAssemblyBegin(*Pv);
833:   VecAssemblyEnd(*Pv);

835:   PetscFree(global_indices);
836:   return(0);
837: }