Actual source code: ispai.c
petsc-3.5.4 2015-05-23
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: */
21: #include <petsc-private/pcimpl.h> /*I "petscpc.h" I*/
22: #include <../src/ksp/pc/impls/spai/petscspai.h>
24: /*
25: These are the SPAI include files
26: */
27: EXTERN_C_BEGIN
28: #define SPAI_USE_MPI /* required for setting SPAI_Comm correctly in basics.h */
29: #include <spai.h>
30: #include <matrix.h>
31: EXTERN_C_END
33: extern PetscErrorCode ConvertMatToMatrix(MPI_Comm,Mat,Mat,matrix**);
34: extern PetscErrorCode ConvertMatrixToMat(MPI_Comm,matrix*,Mat*);
35: extern PetscErrorCode ConvertVectorToVec(MPI_Comm,vector *v,Vec *Pv);
36: extern PetscErrorCode MM_to_PETSC(char*,char*,char*);
38: typedef struct {
40: matrix *B; /* matrix in SPAI format */
41: matrix *BT; /* transpose of matrix in SPAI format */
42: matrix *M; /* the approximate inverse in SPAI format */
44: Mat PM; /* the approximate inverse PETSc format */
46: double epsilon; /* tolerance */
47: int nbsteps; /* max number of "improvement" steps per line */
48: int max; /* max dimensions of is_I, q, etc. */
49: int maxnew; /* max number of new entries per step */
50: int block_size; /* constant block size */
51: int cache_size; /* one of (1,2,3,4,5,6) indicting size of cache */
52: int verbose; /* SPAI prints timing and statistics */
54: int sp; /* symmetric nonzero pattern */
55: MPI_Comm comm_spai; /* communicator to be used with spai */
56: } PC_SPAI;
58: /**********************************************************************/
62: static PetscErrorCode PCSetUp_SPAI(PC pc)
63: {
64: PC_SPAI *ispai = (PC_SPAI*)pc->data;
66: Mat AT;
69: init_SPAI();
71: if (ispai->sp) {
72: ConvertMatToMatrix(ispai->comm_spai,pc->pmat,pc->pmat,&ispai->B);
73: } else {
74: /* Use the transpose to get the column nonzero structure. */
75: MatTranspose(pc->pmat,MAT_INITIAL_MATRIX,&AT);
76: ConvertMatToMatrix(ispai->comm_spai,pc->pmat,AT,&ispai->B);
77: MatDestroy(&AT);
78: }
80: /* Destroy the transpose */
81: /* Don't know how to do it. PETSc developers? */
83: /* construct SPAI preconditioner */
84: /* FILE *messages */ /* file for warning messages */
85: /* double epsilon */ /* tolerance */
86: /* int nbsteps */ /* max number of "improvement" steps per line */
87: /* int max */ /* max dimensions of is_I, q, etc. */
88: /* int maxnew */ /* max number of new entries per step */
89: /* int block_size */ /* block_size == 1 specifies scalar elments
90: block_size == n specifies nxn constant-block elements
91: block_size == 0 specifies variable-block elements */
92: /* int cache_size */ /* one of (1,2,3,4,5,6) indicting size of cache. cache_size == 0 indicates no caching */
93: /* int verbose */ /* verbose == 0 specifies that SPAI is silent
94: verbose == 1 prints timing and matrix statistics */
96: bspai(ispai->B,&ispai->M,
97: stdout,
98: ispai->epsilon,
99: ispai->nbsteps,
100: ispai->max,
101: ispai->maxnew,
102: ispai->block_size,
103: ispai->cache_size,
104: ispai->verbose);
106: ConvertMatrixToMat(PetscObjectComm((PetscObject)pc),ispai->M,&ispai->PM);
108: /* free the SPAI matrices */
109: sp_free_matrix(ispai->B);
110: sp_free_matrix(ispai->M);
111: return(0);
112: }
114: /**********************************************************************/
118: static PetscErrorCode PCApply_SPAI(PC pc,Vec xx,Vec y)
119: {
120: PC_SPAI *ispai = (PC_SPAI*)pc->data;
124: /* Now using PETSc's multiply */
125: MatMult(ispai->PM,xx,y);
126: return(0);
127: }
129: /**********************************************************************/
133: static PetscErrorCode PCDestroy_SPAI(PC pc)
134: {
136: PC_SPAI *ispai = (PC_SPAI*)pc->data;
139: MatDestroy(&ispai->PM);
140: MPI_Comm_free(&(ispai->comm_spai));
141: PetscFree(pc->data);
142: return(0);
143: }
145: /**********************************************************************/
149: static PetscErrorCode PCView_SPAI(PC pc,PetscViewer viewer)
150: {
151: PC_SPAI *ispai = (PC_SPAI*)pc->data;
153: PetscBool iascii;
156: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
157: if (iascii) {
158: PetscViewerASCIIPrintf(viewer," SPAI preconditioner\n");
159: PetscViewerASCIIPrintf(viewer," epsilon %g\n", (double)ispai->epsilon);
160: PetscViewerASCIIPrintf(viewer," nbsteps %d\n", ispai->nbsteps);
161: PetscViewerASCIIPrintf(viewer," max %d\n", ispai->max);
162: PetscViewerASCIIPrintf(viewer," maxnew %d\n", ispai->maxnew);
163: PetscViewerASCIIPrintf(viewer," block_size %d\n",ispai->block_size);
164: PetscViewerASCIIPrintf(viewer," cache_size %d\n",ispai->cache_size);
165: PetscViewerASCIIPrintf(viewer," verbose %d\n", ispai->verbose);
166: PetscViewerASCIIPrintf(viewer," sp %d\n", ispai->sp);
167: }
168: return(0);
169: }
173: static PetscErrorCode PCSPAISetEpsilon_SPAI(PC pc,double epsilon1)
174: {
175: PC_SPAI *ispai = (PC_SPAI*)pc->data;
178: ispai->epsilon = epsilon1;
179: return(0);
180: }
182: /**********************************************************************/
186: static PetscErrorCode PCSPAISetNBSteps_SPAI(PC pc,int nbsteps1)
187: {
188: PC_SPAI *ispai = (PC_SPAI*)pc->data;
191: ispai->nbsteps = nbsteps1;
192: return(0);
193: }
195: /**********************************************************************/
197: /* added 1/7/99 g.h. */
200: static PetscErrorCode PCSPAISetMax_SPAI(PC pc,int max1)
201: {
202: PC_SPAI *ispai = (PC_SPAI*)pc->data;
205: ispai->max = max1;
206: return(0);
207: }
209: /**********************************************************************/
213: static PetscErrorCode PCSPAISetMaxNew_SPAI(PC pc,int maxnew1)
214: {
215: PC_SPAI *ispai = (PC_SPAI*)pc->data;
218: ispai->maxnew = maxnew1;
219: return(0);
220: }
222: /**********************************************************************/
226: static PetscErrorCode PCSPAISetBlockSize_SPAI(PC pc,int block_size1)
227: {
228: PC_SPAI *ispai = (PC_SPAI*)pc->data;
231: ispai->block_size = block_size1;
232: return(0);
233: }
235: /**********************************************************************/
239: static PetscErrorCode PCSPAISetCacheSize_SPAI(PC pc,int cache_size)
240: {
241: PC_SPAI *ispai = (PC_SPAI*)pc->data;
244: ispai->cache_size = cache_size;
245: return(0);
246: }
248: /**********************************************************************/
252: static PetscErrorCode PCSPAISetVerbose_SPAI(PC pc,int verbose)
253: {
254: PC_SPAI *ispai = (PC_SPAI*)pc->data;
257: ispai->verbose = verbose;
258: return(0);
259: }
261: /**********************************************************************/
265: static PetscErrorCode PCSPAISetSp_SPAI(PC pc,int sp)
266: {
267: PC_SPAI *ispai = (PC_SPAI*)pc->data;
270: ispai->sp = sp;
271: return(0);
272: }
274: /* -------------------------------------------------------------------*/
278: /*@
279: PCSPAISetEpsilon -- Set the tolerance for the SPAI preconditioner
281: Input Parameters:
282: + pc - the preconditioner
283: - eps - epsilon (default .4)
285: Notes: Espilon must be between 0 and 1. It controls the
286: quality of the approximation of M to the inverse of
287: A. Higher values of epsilon lead to more work, more
288: fill, and usually better preconditioners. In many
289: cases the best choice of epsilon is the one that
290: divides the total solution time equally between the
291: preconditioner and the solver.
293: Level: intermediate
295: .seealso: PCSPAI, PCSetType()
296: @*/
297: PetscErrorCode PCSPAISetEpsilon(PC pc,double epsilon1)
298: {
302: PetscTryMethod(pc,"PCSPAISetEpsilon_C",(PC,double),(pc,epsilon1));
303: return(0);
304: }
306: /**********************************************************************/
310: /*@
311: PCSPAISetNBSteps - set maximum number of improvement steps per row in
312: the SPAI preconditioner
314: Input Parameters:
315: + pc - the preconditioner
316: - n - number of steps (default 5)
318: Notes: SPAI constructs to approximation to every column of
319: the exact inverse of A in a series of improvement
320: steps. The quality of the approximation is determined
321: by epsilon. If an approximation achieving an accuracy
322: of epsilon is not obtained after ns steps, SPAI simply
323: uses the best approximation constructed so far.
325: Level: intermediate
327: .seealso: PCSPAI, PCSetType(), PCSPAISetMaxNew()
328: @*/
329: PetscErrorCode PCSPAISetNBSteps(PC pc,int nbsteps1)
330: {
334: PetscTryMethod(pc,"PCSPAISetNBSteps_C",(PC,int),(pc,nbsteps1));
335: return(0);
336: }
338: /**********************************************************************/
340: /* added 1/7/99 g.h. */
343: /*@
344: PCSPAISetMax - set the size of various working buffers in
345: the SPAI preconditioner
347: Input Parameters:
348: + pc - the preconditioner
349: - n - size (default is 5000)
351: Level: intermediate
353: .seealso: PCSPAI, PCSetType()
354: @*/
355: PetscErrorCode PCSPAISetMax(PC pc,int max1)
356: {
360: PetscTryMethod(pc,"PCSPAISetMax_C",(PC,int),(pc,max1));
361: return(0);
362: }
364: /**********************************************************************/
368: /*@
369: PCSPAISetMaxNew - set maximum number of new nonzero candidates per step
370: in SPAI preconditioner
372: Input Parameters:
373: + pc - the preconditioner
374: - n - maximum number (default 5)
376: Level: intermediate
378: .seealso: PCSPAI, PCSetType(), PCSPAISetNBSteps()
379: @*/
380: PetscErrorCode PCSPAISetMaxNew(PC pc,int maxnew1)
381: {
385: PetscTryMethod(pc,"PCSPAISetMaxNew_C",(PC,int),(pc,maxnew1));
386: return(0);
387: }
389: /**********************************************************************/
393: /*@
394: PCSPAISetBlockSize - set the block size for the SPAI preconditioner
396: Input Parameters:
397: + pc - the preconditioner
398: - n - block size (default 1)
400: Notes: A block
401: size of 1 treats A as a matrix of scalar elements. A
402: block size of s > 1 treats A as a matrix of sxs
403: blocks. A block size of 0 treats A as a matrix with
404: variable sized blocks, which are determined by
405: searching for dense square diagonal blocks in A.
406: This can be very effective for finite-element
407: matrices.
409: SPAI will convert A to block form, use a block
410: version of the preconditioner algorithm, and then
411: convert the result back to scalar form.
413: In many cases the a block-size parameter other than 1
414: can lead to very significant improvement in
415: performance.
418: Level: intermediate
420: .seealso: PCSPAI, PCSetType()
421: @*/
422: PetscErrorCode PCSPAISetBlockSize(PC pc,int block_size1)
423: {
427: PetscTryMethod(pc,"PCSPAISetBlockSize_C",(PC,int),(pc,block_size1));
428: return(0);
429: }
431: /**********************************************************************/
435: /*@
436: PCSPAISetCacheSize - specify cache size in the SPAI preconditioner
438: Input Parameters:
439: + pc - the preconditioner
440: - n - cache size {0,1,2,3,4,5} (default 5)
442: Notes: SPAI uses a hash table to cache messages and avoid
443: redundant communication. If suggest always using
444: 5. This parameter is irrelevant in the serial
445: version.
447: Level: intermediate
449: .seealso: PCSPAI, PCSetType()
450: @*/
451: PetscErrorCode PCSPAISetCacheSize(PC pc,int cache_size)
452: {
456: PetscTryMethod(pc,"PCSPAISetCacheSize_C",(PC,int),(pc,cache_size));
457: return(0);
458: }
460: /**********************************************************************/
464: /*@
465: PCSPAISetVerbose - verbosity level for the SPAI preconditioner
467: Input Parameters:
468: + pc - the preconditioner
469: - n - level (default 1)
471: Notes: print parameters, timings and matrix statistics
473: Level: intermediate
475: .seealso: PCSPAI, PCSetType()
476: @*/
477: PetscErrorCode PCSPAISetVerbose(PC pc,int verbose)
478: {
482: PetscTryMethod(pc,"PCSPAISetVerbose_C",(PC,int),(pc,verbose));
483: return(0);
484: }
486: /**********************************************************************/
490: /*@
491: PCSPAISetSp - specify a symmetric matrix sparsity pattern in the SPAI preconditioner
493: Input Parameters:
494: + pc - the preconditioner
495: - n - 0 or 1
497: Notes: If A has a symmetric nonzero pattern use -sp 1 to
498: improve performance by eliminating some communication
499: in the parallel version. Even if A does not have a
500: symmetric nonzero pattern -sp 1 may well lead to good
501: results, but the code will not follow the published
502: SPAI algorithm exactly.
505: Level: intermediate
507: .seealso: PCSPAI, PCSetType()
508: @*/
509: PetscErrorCode PCSPAISetSp(PC pc,int sp)
510: {
514: PetscTryMethod(pc,"PCSPAISetSp_C",(PC,int),(pc,sp));
515: return(0);
516: }
518: /**********************************************************************/
520: /**********************************************************************/
524: static PetscErrorCode PCSetFromOptions_SPAI(PC pc)
525: {
526: PC_SPAI *ispai = (PC_SPAI*)pc->data;
528: int nbsteps1,max1,maxnew1,block_size1,cache_size,verbose,sp;
529: double epsilon1;
530: PetscBool flg;
533: PetscOptionsHead("SPAI options");
534: PetscOptionsReal("-pc_spai_epsilon","","PCSPAISetEpsilon",ispai->epsilon,&epsilon1,&flg);
535: if (flg) {
536: PCSPAISetEpsilon(pc,epsilon1);
537: }
538: PetscOptionsInt("-pc_spai_nbsteps","","PCSPAISetNBSteps",ispai->nbsteps,&nbsteps1,&flg);
539: if (flg) {
540: PCSPAISetNBSteps(pc,nbsteps1);
541: }
542: /* added 1/7/99 g.h. */
543: PetscOptionsInt("-pc_spai_max","","PCSPAISetMax",ispai->max,&max1,&flg);
544: if (flg) {
545: PCSPAISetMax(pc,max1);
546: }
547: PetscOptionsInt("-pc_spai_maxnew","","PCSPAISetMaxNew",ispai->maxnew,&maxnew1,&flg);
548: if (flg) {
549: PCSPAISetMaxNew(pc,maxnew1);
550: }
551: PetscOptionsInt("-pc_spai_block_size","","PCSPAISetBlockSize",ispai->block_size,&block_size1,&flg);
552: if (flg) {
553: PCSPAISetBlockSize(pc,block_size1);
554: }
555: PetscOptionsInt("-pc_spai_cache_size","","PCSPAISetCacheSize",ispai->cache_size,&cache_size,&flg);
556: if (flg) {
557: PCSPAISetCacheSize(pc,cache_size);
558: }
559: PetscOptionsInt("-pc_spai_verbose","","PCSPAISetVerbose",ispai->verbose,&verbose,&flg);
560: if (flg) {
561: PCSPAISetVerbose(pc,verbose);
562: }
563: PetscOptionsInt("-pc_spai_sp","","PCSPAISetSp",ispai->sp,&sp,&flg);
564: if (flg) {
565: PCSPAISetSp(pc,sp);
566: }
567: PetscOptionsTail();
568: return(0);
569: }
571: /**********************************************************************/
573: /*MC
574: PCSPAI - Use the Sparse Approximate Inverse method of Grote and Barnard
575: as a preconditioner (SIAM J. Sci. Comput.; vol 18, nr 3)
577: Options Database Keys:
578: + -pc_spai_epsilon <eps> - set tolerance
579: . -pc_spai_nbstep <n> - set nbsteps
580: . -pc_spai_max <m> - set max
581: . -pc_spai_max_new <m> - set maxnew
582: . -pc_spai_block_size <n> - set block size
583: . -pc_spai_cache_size <n> - set cache size
584: . -pc_spai_sp <m> - set sp
585: - -pc_spai_set_verbose <true,false> - verbose output
587: Notes: This only works with AIJ matrices.
589: Level: beginner
591: Concepts: approximate inverse
593: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
594: PCSPAISetEpsilon(), PCSPAISetMax(), PCSPAISetMaxNew(), PCSPAISetBlockSize(),
595: PCSPAISetVerbose(), PCSPAISetSp()
596: M*/
600: PETSC_EXTERN PetscErrorCode PCCreate_SPAI(PC pc)
601: {
602: PC_SPAI *ispai;
606: PetscNewLog(pc,&ispai);
607: pc->data = ispai;
609: pc->ops->destroy = PCDestroy_SPAI;
610: pc->ops->apply = PCApply_SPAI;
611: pc->ops->applyrichardson = 0;
612: pc->ops->setup = PCSetUp_SPAI;
613: pc->ops->view = PCView_SPAI;
614: pc->ops->setfromoptions = PCSetFromOptions_SPAI;
616: ispai->epsilon = .4;
617: ispai->nbsteps = 5;
618: ispai->max = 5000;
619: ispai->maxnew = 5;
620: ispai->block_size = 1;
621: ispai->cache_size = 5;
622: ispai->verbose = 0;
624: ispai->sp = 1;
625: MPI_Comm_dup(PetscObjectComm((PetscObject)pc),&(ispai->comm_spai));
627: PetscObjectComposeFunction((PetscObject)pc,"PCSPAISetEpsilon_C",PCSPAISetEpsilon_SPAI);
628: PetscObjectComposeFunction((PetscObject)pc,"PCSPAISetNBSteps_C",PCSPAISetNBSteps_SPAI);
629: PetscObjectComposeFunction((PetscObject)pc,"PCSPAISetMax_C",PCSPAISetMax_SPAI);
630: PetscObjectComposeFunction((PetscObject)pc,"PCSPAISetMaxNew_C",PCSPAISetMaxNew_SPAI);
631: PetscObjectComposeFunction((PetscObject)pc,"PCSPAISetBlockSize_C",PCSPAISetBlockSize_SPAI);
632: PetscObjectComposeFunction((PetscObject)pc,"PCSPAISetCacheSize_C",PCSPAISetCacheSize_SPAI);
633: PetscObjectComposeFunction((PetscObject)pc,"PCSPAISetVerbose_C",PCSPAISetVerbose_SPAI);
634: PetscObjectComposeFunction((PetscObject)pc,"PCSPAISetSp_C",PCSPAISetSp_SPAI);
635: return(0);
636: }
638: /**********************************************************************/
640: /*
641: Converts from a PETSc matrix to an SPAI matrix
642: */
645: PetscErrorCode ConvertMatToMatrix(MPI_Comm comm, Mat A,Mat AT,matrix **B)
646: {
647: matrix *M;
648: int i,j,col;
649: int row_indx;
650: int len,pe,local_indx,start_indx;
651: int *mapping;
652: PetscErrorCode ierr;
653: const int *cols;
654: const double *vals;
655: int n,mnl,nnl,nz,rstart,rend;
656: PetscMPIInt size,rank;
657: struct compressed_lines *rows;
660: MPI_Comm_size(comm,&size);
661: MPI_Comm_rank(comm,&rank);
662: MatGetSize(A,&n,&n);
663: MatGetLocalSize(A,&mnl,&nnl);
665: /*
666: not sure why a barrier is required. commenting out
667: MPI_Barrier(comm);
668: */
670: M = new_matrix((SPAI_Comm)comm);
672: M->n = n;
673: M->bs = 1;
674: M->max_block_size = 1;
676: M->mnls = (int*)malloc(sizeof(int)*size);
677: M->start_indices = (int*)malloc(sizeof(int)*size);
678: M->pe = (int*)malloc(sizeof(int)*n);
679: M->block_sizes = (int*)malloc(sizeof(int)*n);
680: for (i=0; i<n; i++) M->block_sizes[i] = 1;
682: MPI_Allgather(&mnl,1,MPI_INT,M->mnls,1,MPI_INT,comm);
684: M->start_indices[0] = 0;
685: for (i=1; i<size; i++) M->start_indices[i] = M->start_indices[i-1] + M->mnls[i-1];
687: M->mnl = M->mnls[M->myid];
688: M->my_start_index = M->start_indices[M->myid];
690: for (i=0; i<size; i++) {
691: start_indx = M->start_indices[i];
692: for (j=0; j<M->mnls[i]; j++) M->pe[start_indx+j] = i;
693: }
695: if (AT) {
696: M->lines = new_compressed_lines(M->mnls[rank],1);
697: } else {
698: M->lines = new_compressed_lines(M->mnls[rank],0);
699: }
701: rows = M->lines;
703: /* Determine the mapping from global indices to pointers */
704: PetscMalloc1(M->n,&mapping);
705: pe = 0;
706: local_indx = 0;
707: for (i=0; i<M->n; i++) {
708: if (local_indx >= M->mnls[pe]) {
709: pe++;
710: local_indx = 0;
711: }
712: mapping[i] = local_indx + M->start_indices[pe];
713: local_indx++;
714: }
716: /*********************************************************/
717: /************** Set up the row structure *****************/
718: /*********************************************************/
720: MatGetOwnershipRange(A,&rstart,&rend);
721: for (i=rstart; i<rend; i++) {
722: row_indx = i - rstart;
723: MatGetRow(A,i,&nz,&cols,&vals);
724: /* allocate buffers */
725: rows->ptrs[row_indx] = (int*)malloc(nz*sizeof(int));
726: rows->A[row_indx] = (double*)malloc(nz*sizeof(double));
727: /* copy the matrix */
728: for (j=0; j<nz; j++) {
729: col = cols[j];
730: len = rows->len[row_indx]++;
732: rows->ptrs[row_indx][len] = mapping[col];
733: rows->A[row_indx][len] = vals[j];
734: }
735: rows->slen[row_indx] = rows->len[row_indx];
737: MatRestoreRow(A,i,&nz,&cols,&vals);
738: }
741: /************************************************************/
742: /************** Set up the column structure *****************/
743: /*********************************************************/
745: if (AT) {
747: for (i=rstart; i<rend; i++) {
748: row_indx = i - rstart;
749: MatGetRow(AT,i,&nz,&cols,&vals);
750: /* allocate buffers */
751: rows->rptrs[row_indx] = (int*)malloc(nz*sizeof(int));
752: /* copy the matrix (i.e., the structure) */
753: for (j=0; j<nz; j++) {
754: col = cols[j];
755: len = rows->rlen[row_indx]++;
757: rows->rptrs[row_indx][len] = mapping[col];
758: }
759: MatRestoreRow(AT,i,&nz,&cols,&vals);
760: }
761: }
763: PetscFree(mapping);
765: order_pointers(M);
766: M->maxnz = calc_maxnz(M);
767: *B = M;
768: return(0);
769: }
771: /**********************************************************************/
773: /*
774: Converts from an SPAI matrix B to a PETSc matrix PB.
775: This assumes that the SPAI matrix B is stored in
776: COMPRESSED-ROW format.
777: */
780: PetscErrorCode ConvertMatrixToMat(MPI_Comm comm,matrix *B,Mat *PB)
781: {
782: PetscMPIInt size,rank;
784: int m,n,M,N;
785: int d_nz,o_nz;
786: int *d_nnz,*o_nnz;
787: int i,k,global_row,global_col,first_diag_col,last_diag_col;
788: PetscScalar val;
791: MPI_Comm_size(comm,&size);
792: MPI_Comm_rank(comm,&rank);
794: m = n = B->mnls[rank];
795: d_nz = o_nz = 0;
797: /* Determine preallocation for MatCreateMPIAIJ */
798: PetscMalloc1(m,&d_nnz);
799: PetscMalloc1(m,&o_nnz);
800: for (i=0; i<m; i++) d_nnz[i] = o_nnz[i] = 0;
801: first_diag_col = B->start_indices[rank];
802: last_diag_col = first_diag_col + B->mnls[rank];
803: for (i=0; i<B->mnls[rank]; i++) {
804: for (k=0; k<B->lines->len[i]; k++) {
805: global_col = B->lines->ptrs[i][k];
806: if ((global_col >= first_diag_col) && (global_col < last_diag_col)) d_nnz[i]++;
807: else o_nnz[i]++;
808: }
809: }
811: M = N = B->n;
812: /* Here we only know how to create AIJ format */
813: MatCreate(comm,PB);
814: MatSetSizes(*PB,m,n,M,N);
815: MatSetType(*PB,MATAIJ);
816: MatSeqAIJSetPreallocation(*PB,d_nz,d_nnz);
817: MatMPIAIJSetPreallocation(*PB,d_nz,d_nnz,o_nz,o_nnz);
819: for (i=0; i<B->mnls[rank]; i++) {
820: global_row = B->start_indices[rank]+i;
821: for (k=0; k<B->lines->len[i]; k++) {
822: global_col = B->lines->ptrs[i][k];
824: val = B->lines->A[i][k];
825: MatSetValues(*PB,1,&global_row,1,&global_col,&val,ADD_VALUES);
826: }
827: }
829: PetscFree(d_nnz);
830: PetscFree(o_nnz);
832: MatAssemblyBegin(*PB,MAT_FINAL_ASSEMBLY);
833: MatAssemblyEnd(*PB,MAT_FINAL_ASSEMBLY);
834: return(0);
835: }
837: /**********************************************************************/
839: /*
840: Converts from an SPAI vector v to a PETSc vec Pv.
841: */
844: PetscErrorCode ConvertVectorToVec(MPI_Comm comm,vector *v,Vec *Pv)
845: {
847: PetscMPIInt size,rank;
848: int m,M,i,*mnls,*start_indices,*global_indices;
851: MPI_Comm_size(comm,&size);
852: MPI_Comm_rank(comm,&rank);
854: m = v->mnl;
855: M = v->n;
858: VecCreateMPI(comm,m,M,Pv);
860: PetscMalloc1(size,&mnls);
861: MPI_Allgather(&v->mnl,1,MPI_INT,mnls,1,MPI_INT,comm);
863: PetscMalloc1(size,&start_indices);
865: start_indices[0] = 0;
866: for (i=1; i<size; i++) start_indices[i] = start_indices[i-1] +mnls[i-1];
868: PetscMalloc1(v->mnl,&global_indices);
869: for (i=0; i<v->mnl; i++) global_indices[i] = start_indices[rank] + i;
871: PetscFree(mnls);
872: PetscFree(start_indices);
874: VecSetValues(*Pv,v->mnl,global_indices,v->v,INSERT_VALUES);
875: VecAssemblyBegin(*Pv);
876: VecAssemblyEnd(*Pv);
878: PetscFree(global_indices);
879: return(0);
880: }