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: }