Actual source code: plapack.c
1: #define PETSCMAT_DLL
3: /*
4: Provides an interface to the PLAPACKR32 dense solver
5: */
7: #include src/mat/impls/dense/seq/dense.h
8: #include src/mat/impls/dense/mpi/mpidense.h
11: #include "PLA.h"
12: #include "PLA_prototypes.h"
15: typedef struct {
16: MPI_Comm comm_2d;
17: PLA_Obj A,pivots;
18: PLA_Template templ;
19: MPI_Datatype datatype;
20: PetscInt nb,nb_alg,ierror,rstart;
21: VecScatter ctx;
22: IS is_pla,is_petsc;
23: PetscTruth pla_solved;
24: MatStructure mstruct;
25: PetscMPIInt nprows,npcols;
27: /* A few function pointers for inheritance */
28: PetscErrorCode (*MatDuplicate)(Mat,MatDuplicateOption,Mat*);
29: PetscErrorCode (*MatView)(Mat,PetscViewer);
30: PetscErrorCode (*MatAssemblyEnd)(Mat,MatAssemblyType);
31: PetscErrorCode (*MatLUFactorSymbolic)(Mat,IS,IS,MatFactorInfo*,Mat*);
32: PetscErrorCode (*MatCholeskyFactorSymbolic)(Mat,IS,MatFactorInfo*,Mat*);
33: PetscErrorCode (*MatDestroy)(Mat);
35: /* Flag to clean up (non-global) Plapack objects during Destroy */
36: PetscTruth CleanUpPlapack;
37: } Mat_Plapack;
39: EXTERN PetscErrorCode MatDuplicate_Plapack(Mat,MatDuplicateOption,Mat*);
44: PetscErrorCode MatConvert_Plapack_Dense(Mat A,MatType type,MatReuse reuse,Mat *newmat)
45: {
47: PetscErrorCode ierr;
48: Mat B=*newmat;
49: Mat_Plapack *lu=(Mat_Plapack *)A->spptr;
52: if (reuse == MAT_INITIAL_MATRIX) {
53: MatDuplicate(A,MAT_COPY_VALUES,&B);
54: }
55: /* Reset the original function pointers */
56: B->ops->duplicate = lu->MatDuplicate;
57: B->ops->view = lu->MatView;
58: B->ops->assemblyend = lu->MatAssemblyEnd;
59: B->ops->lufactorsymbolic = lu->MatLUFactorSymbolic;
60: B->ops->destroy = lu->MatDestroy;
62: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_plapack_C","",PETSC_NULL);
63: PetscObjectComposeFunction((PetscObject)B,"MatConvert_plapack_seqdense_C","",PETSC_NULL);
64: PetscObjectComposeFunction((PetscObject)B,"MatConvert_mpidense_plapack_C","",PETSC_NULL);
65: PetscObjectComposeFunction((PetscObject)B,"MatConvert_plapack_mpidense_C","",PETSC_NULL);
67: PetscObjectChangeTypeName((PetscObject)B,type);
68: *newmat = B;
69: return(0);
70: }
75: PetscErrorCode MatDestroy_Plapack(Mat A)
76: {
78: PetscMPIInt size;
79: Mat_Plapack *lu=(Mat_Plapack*)A->spptr;
80:
82: if (lu->CleanUpPlapack) {
83: /* Deallocate Plapack storage */
84: PLA_Obj_free(&lu->A);
85: PLA_Obj_free (&lu->pivots);
86: PLA_Temp_free(&lu->templ);
87: PLA_Finalize();
89: ISDestroy(lu->is_pla);
90: ISDestroy(lu->is_petsc);
91: VecScatterDestroy(lu->ctx);
92: }
93: MPI_Comm_size(A->comm,&size);
94: if (size == 1) {
95: MatConvert_Plapack_Dense(A,MATSEQDENSE,MAT_REUSE_MATRIX,&A);
96: } else {
97: MatConvert_Plapack_Dense(A,MATMPIDENSE,MAT_REUSE_MATRIX,&A);
98: }
99: (*A->ops->destroy)(A);
100: return(0);
101: }
105: PetscErrorCode MatSolve_Plapack(Mat A,Vec b,Vec x)
106: {
107: MPI_Comm comm = A->comm;
108: Mat_Plapack *lu = (Mat_Plapack*)A->spptr;
110: PetscInt M=A->rmap.N,m=A->rmap.n,rstart,i,j,*idx_pla,*idx_petsc,loc_m,loc_stride;
111: PetscScalar *array;
112: PetscReal one = 1.0;
113: PetscMPIInt size,rank,r_rank,r_nproc,c_rank,c_nproc;;
114: PLA_Obj v_pla = NULL;
115: PetscScalar *loc_buf;
116: Vec loc_x;
117:
119: MPI_Comm_size(comm,&size);
120: MPI_Comm_rank(comm,&rank);
122: /* Create PLAPACK vector objects, then copy b into PLAPACK b */
123: PLA_Mvector_create(lu->datatype,M,1,lu->templ,PLA_ALIGN_FIRST,&v_pla);
124: PLA_Obj_set_to_zero(v_pla);
126: /* Copy b into rhs_pla */
127: PLA_API_begin();
128: PLA_Obj_API_open(v_pla);
129: VecGetArray(b,&array);
130: PLA_API_axpy_vector_to_global(m,&one,(void *)array,1,v_pla,lu->rstart);
131: VecRestoreArray(b,&array);
132: PLA_Obj_API_close(v_pla);
133: PLA_API_end();
135: if (A->factor == FACTOR_LU){
136: /* Apply the permutations to the right hand sides */
137: PLA_Apply_pivots_to_rows (v_pla,lu->pivots);
139: /* Solve L y = b, overwriting b with y */
140: PLA_Trsv( PLA_LOWER_TRIANGULAR,PLA_NO_TRANSPOSE,PLA_UNIT_DIAG,lu->A,v_pla );
142: /* Solve U x = y (=b), overwriting b with x */
143: PLA_Trsv( PLA_UPPER_TRIANGULAR,PLA_NO_TRANSPOSE,PLA_NONUNIT_DIAG,lu->A,v_pla );
144: } else { /* FACTOR_CHOLESKY */
145: PLA_Trsv( PLA_LOWER_TRIANGULAR,PLA_NO_TRANSPOSE,PLA_NONUNIT_DIAG,lu->A,v_pla);
146: PLA_Trsv( PLA_LOWER_TRIANGULAR,(lu->datatype == MPI_DOUBLE ? PLA_TRANSPOSE : PLA_CONJUGATE_TRANSPOSE),
147: PLA_NONUNIT_DIAG,lu->A,v_pla);
148: }
150: /* Copy PLAPACK x into Petsc vector x */
151: PLA_Obj_local_length(v_pla, &loc_m);
152: PLA_Obj_local_buffer(v_pla, (void**)&loc_buf);
153: PLA_Obj_local_stride(v_pla, &loc_stride);
154: /*
155: PetscPrintf(PETSC_COMM_SELF," [%d] b - local_m %d local_stride %d, loc_buf: %g %g, nb: %d\n",rank,loc_m,loc_stride,loc_buf[0],loc_buf[(loc_m-1)*loc_stride],lu->nb);
156: */
157: VecCreateSeqWithArray(PETSC_COMM_SELF,loc_m*loc_stride,loc_buf,&loc_x);
158: if (!lu->pla_solved){
159:
160: PLA_Temp_comm_row_info(lu->templ,&lu->comm_2d,&r_rank,&r_nproc);
161: PLA_Temp_comm_col_info(lu->templ,&lu->comm_2d,&c_rank,&c_nproc);
162: /* printf(" [%d] rank: %d %d, nproc: %d %d\n",rank,r_rank,c_rank,r_nproc,c_nproc); */
164: /* Create IS and cts for VecScatterring */
165: PLA_Obj_local_length(v_pla, &loc_m);
166: PLA_Obj_local_stride(v_pla, &loc_stride);
167: PetscMalloc((2*loc_m+1)*sizeof(PetscInt),&idx_pla);
168: idx_petsc = idx_pla + loc_m;
170: rstart = (r_rank*c_nproc+c_rank)*lu->nb;
171: for (i=0; i<loc_m; i+=lu->nb){
172: j = 0;
173: while (j < lu->nb && i+j < loc_m){
174: idx_petsc[i+j] = rstart + j; j++;
175: }
176: rstart += size*lu->nb;
177: }
179: for (i=0; i<loc_m; i++) idx_pla[i] = i*loc_stride;
181: ISCreateGeneral(PETSC_COMM_SELF,loc_m,idx_pla,&lu->is_pla);
182: ISCreateGeneral(PETSC_COMM_SELF,loc_m,idx_petsc,&lu->is_petsc);
183: PetscFree(idx_pla);
184: VecScatterCreate(loc_x,lu->is_pla,x,lu->is_petsc,&lu->ctx);
185: }
186: VecScatterBegin(loc_x,x,INSERT_VALUES,SCATTER_FORWARD,lu->ctx);
187: VecScatterEnd(loc_x,x,INSERT_VALUES,SCATTER_FORWARD,lu->ctx);
188:
189: /* Free data */
190: VecDestroy(loc_x);
191: PLA_Obj_free(&v_pla);
193: lu->pla_solved = PETSC_TRUE;
194: return(0);
195: }
199: PetscErrorCode MatLUFactorNumeric_Plapack(Mat A,MatFactorInfo *info,Mat *F)
200: {
201: Mat_Plapack *lu = (Mat_Plapack*)(*F)->spptr;
203: PetscInt M=A->rmap.N,m=A->rmap.n,rstart,rend;
204: PetscInt info_pla=0;
205: PetscScalar *array,one = 1.0;
208: if (lu->mstruct == SAME_NONZERO_PATTERN){
209: PLA_Obj_free(&lu->A);
210: PLA_Obj_free (&lu->pivots);
211: }
212: /* Create PLAPACK matrix object */
213: lu->A = NULL; lu->pivots = NULL;
214: PLA_Matrix_create(lu->datatype,M,M,lu->templ,PLA_ALIGN_FIRST,PLA_ALIGN_FIRST,&lu->A);
215: PLA_Obj_set_to_zero(lu->A);
216: PLA_Mvector_create(MPI_INT,M,1,lu->templ,PLA_ALIGN_FIRST,&lu->pivots);
218: /* Copy A into lu->A */
219: PLA_API_begin();
220: PLA_Obj_API_open(lu->A);
221: MatGetOwnershipRange(A,&rstart,&rend);
222: MatGetArray(A,&array);
223: PLA_API_axpy_matrix_to_global(m,M, &one,(void *)array,m,lu->A,rstart,0);
224: MatRestoreArray(A,&array);
225: PLA_Obj_API_close(lu->A);
226: PLA_API_end();
228: /* Factor P A -> L U overwriting lower triangular portion of A with L, upper, U */
229: info_pla = PLA_LU(lu->A,lu->pivots);
230: if (info_pla != 0)
231: SETERRQ1(PETSC_ERR_MAT_LU_ZRPVT,"Zero pivot encountered at row %d from PLA_LU()",info_pla);
233: lu->CleanUpPlapack = PETSC_TRUE;
234: lu->rstart = rstart;
235: lu->mstruct = SAME_NONZERO_PATTERN;
236:
237: (*F)->assembled = PETSC_TRUE; /* required by -ksp_view */
238: return(0);
239: }
243: PetscErrorCode MatCholeskyFactorNumeric_Plapack(Mat A,MatFactorInfo *info,Mat *F)
244: {
245: Mat_Plapack *lu = (Mat_Plapack*)(*F)->spptr;
247: PetscInt M=A->rmap.N,m=A->rmap.n,rstart,rend;
248: PetscInt info_pla=0;
249: PetscScalar *array,one = 1.0;
252: if (lu->mstruct == SAME_NONZERO_PATTERN){
253: PLA_Obj_free(&lu->A);
254: }
255: /* Create PLAPACK matrix object */
256: lu->A = NULL;
257: lu->pivots = NULL;
258: PLA_Matrix_create(lu->datatype,M,M,lu->templ,PLA_ALIGN_FIRST,PLA_ALIGN_FIRST,&lu->A);
260: /* Copy A into lu->A */
261: PLA_API_begin();
262: PLA_Obj_API_open(lu->A);
263: MatGetOwnershipRange(A,&rstart,&rend);
264: MatGetArray(A,&array);
265: PLA_API_axpy_matrix_to_global(m,M, &one,(void *)array,m,lu->A,rstart,0);
266: MatRestoreArray(A,&array);
267: PLA_Obj_API_close(lu->A);
268: PLA_API_end();
270: /* Factor P A -> Chol */
271: info_pla = PLA_Chol(PLA_LOWER_TRIANGULAR,lu->A);
272: if (info_pla != 0)
273: SETERRQ1( PETSC_ERR_MAT_CH_ZRPVT,"Nonpositive definite matrix detected at row %d from PLA_Chol()",info_pla);
275: lu->CleanUpPlapack = PETSC_TRUE;
276: lu->rstart = rstart;
277: lu->mstruct = SAME_NONZERO_PATTERN;
278:
279: (*F)->assembled = PETSC_TRUE; /* required by -ksp_view */
280: return(0);
281: }
285: PetscErrorCode MatFactorSymbolic_Plapack_Private(Mat A,MatFactorInfo *info,Mat *F)
286: {
287: Mat B;
288: Mat_Plapack *lu;
290: PetscInt M=A->rmap.N,N=A->cmap.N;
291: MPI_Comm comm=A->comm,comm_2d;
292: PetscMPIInt size;
293: PetscInt ierror;
296: /* Create the factorization matrix */
297: MatCreate(A->comm,&B);
298: MatSetSizes(B,A->rmap.n,A->cmap.n,M,N);
299: MatSetType(B,A->type_name);
301: B->ops->solve = MatSolve_Plapack;
302: lu = (Mat_Plapack*)(B->spptr);
304: /* Set default Plapack parameters */
305: MPI_Comm_size(comm,&size);
306: lu->nprows = 1; lu->npcols = size;
307: ierror = 0;
308: lu->nb = M/size;
309: if (M - lu->nb*size) lu->nb++; /* without cyclic distribution */
310:
311: /* Set runtime options */
312: PetscOptionsBegin(A->comm,A->prefix,"PLAPACK Options","Mat");
313: PetscOptionsInt("-mat_plapack_nprows","row dimension of 2D processor mesh","None",lu->nprows,&lu->nprows,PETSC_NULL);
314: PetscOptionsInt("-mat_plapack_npcols","column dimension of 2D processor mesh","None",lu->npcols,&lu->npcols,PETSC_NULL);
315:
316: PetscOptionsInt("-mat_plapack_nb","block size of template vector","None",lu->nb,&lu->nb,PETSC_NULL);
317: PetscOptionsInt("-mat_plapack_ckerror","error checking flag","None",ierror,&ierror,PETSC_NULL);
318: if (ierror){
319: PLA_Set_error_checking(ierror,PETSC_TRUE,PETSC_TRUE,PETSC_FALSE );
320: } else {
321: PLA_Set_error_checking(ierror,PETSC_FALSE,PETSC_FALSE,PETSC_FALSE );
322: }
323: lu->ierror = ierror;
324:
325: lu->nb_alg = 0;
326: PetscOptionsInt("-mat_plapack_nb_alg","algorithmic block size","None",lu->nb_alg,&lu->nb_alg,PETSC_NULL);
327: if (lu->nb_alg){
328: pla_Environ_set_nb_alg (PLA_OP_ALL_ALG,lu->nb_alg);
329: }
330: PetscOptionsEnd();
333: /* Create a 2D communicator */
334: PLA_Comm_1D_to_2D(comm,lu->nprows,lu->npcols,&comm_2d);
335: lu->comm_2d = comm_2d;
337: /* Initialize PLAPACK */
338: PLA_Init(comm_2d);
340: /* Create object distribution template */
341: lu->templ = NULL;
342: PLA_Temp_create(lu->nb, 0, &lu->templ);
344: /* Use suggested nb_alg if it is not provided by user */
345: if (lu->nb_alg == 0){
346: PLA_Environ_nb_alg(PLA_OP_PAN_PAN,lu->templ,&lu->nb_alg);
347: pla_Environ_set_nb_alg(PLA_OP_ALL_ALG,lu->nb_alg);
348: }
350: /* Set the datatype */
351: #if defined(PETSC_USE_COMPLEX)
352: lu->datatype = MPI_DOUBLE_COMPLEX;
353: #else
354: lu->datatype = MPI_DOUBLE;
355: #endif
357: lu->pla_solved = PETSC_FALSE; /* MatSolve_Plapack() is called yet */
358: lu->mstruct = DIFFERENT_NONZERO_PATTERN;
359: lu->CleanUpPlapack = PETSC_TRUE;
360: *F = B;
361: return(0);
362: }
364: /* Note the Petsc r and c permutations are ignored */
367: PetscErrorCode MatLUFactorSymbolic_Plapack(Mat A,IS r,IS c,MatFactorInfo *info,Mat *F)
368: {
372: MatFactorSymbolic_Plapack_Private(A,info,F);
373: (*F)->ops->lufactornumeric = MatLUFactorNumeric_Plapack;
374: (*F)->factor = FACTOR_LU;
375: return(0);
376: }
378: /* Note the Petsc perm permutation is ignored */
381: PetscErrorCode MatCholeskyFactorSymbolic_Plapack(Mat A,IS perm,MatFactorInfo *info,Mat *F)
382: {
384: PetscTruth issymmetric,set;
387: MatIsSymmetricKnown(A,&set,&issymmetric);
388: if (!set || !issymmetric) SETERRQ(PETSC_ERR_USER,"Matrix must be set as MAT_SYMMETRIC for CholeskyFactor()");
389: MatFactorSymbolic_Plapack_Private(A,info,F);
390: (*F)->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_Plapack;
391: (*F)->factor = FACTOR_CHOLESKY;
392: return(0);
393: }
397: PetscErrorCode MatAssemblyEnd_Plapack(Mat A,MatAssemblyType mode)
398: {
399: PetscErrorCode ierr;
400: Mat_Plapack *lu=(Mat_Plapack*)(A->spptr);
403: (*lu->MatAssemblyEnd)(A,mode);
404: lu->MatLUFactorSymbolic = A->ops->lufactorsymbolic;
405: A->ops->lufactorsymbolic = MatLUFactorSymbolic_Plapack;
406: lu->MatCholeskyFactorSymbolic = A->ops->choleskyfactorsymbolic;
407: A->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_Plapack;
408: return(0);
409: }
413: PetscErrorCode MatFactorInfo_Plapack(Mat A,PetscViewer viewer)
414: {
415: Mat_Plapack *lu=(Mat_Plapack*)A->spptr;
416: PetscErrorCode ierr;
419: /* check if matrix is plapack type */
420: if (A->ops->solve != MatSolve_Plapack) return(0);
422: PetscViewerASCIIPrintf(viewer,"PLAPACK run parameters:\n");
423: PetscViewerASCIIPrintf(viewer," Processor mesh: nprows %d, npcols %d\n",lu->nprows, lu->npcols);
424: PetscViewerASCIIPrintf(viewer," Distr. block size nb: %d \n",lu->nb);
425: PetscViewerASCIIPrintf(viewer," Error checking: %d\n",lu->ierror);
426: PetscViewerASCIIPrintf(viewer," Algorithmic block size: %d\n",lu->nb_alg);
427: return(0);
428: }
432: PetscErrorCode MatView_Plapack(Mat A,PetscViewer viewer)
433: {
434: PetscErrorCode ierr;
435: PetscTruth iascii;
436: PetscViewerFormat format;
437: /* Mat_Plapack *lu=(Mat_Plapack*)(A->spptr); */
440: /* (*lu->MatView)(A,viewer); MatView_MPIDense() crash! */
441: PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);
442: if (iascii) {
443: PetscViewerGetFormat(viewer,&format);
444: if (format == PETSC_VIEWER_ASCII_INFO) {
445: MatFactorInfo_Plapack(A,viewer);
446: }
447: }
448: return(0);
449: }
454: PetscErrorCode MatConvert_Dense_Plapack(Mat A,MatType type,MatReuse reuse,Mat *newmat)
455: {
456: /* This routine is only called to convert to MATPLAPACK from MATDENSE, so we ignore 'MatType type'. */
458: PetscMPIInt size;
459: Mat B=*newmat;
460: Mat_Plapack *lu;
463: if (reuse == MAT_INITIAL_MATRIX) {
464: MatDuplicate(A,MAT_COPY_VALUES,&B);
465: }
467: PetscNew(Mat_Plapack,&lu);
468: lu->MatDuplicate = A->ops->duplicate;
469: lu->MatView = A->ops->view;
470: lu->MatAssemblyEnd = A->ops->assemblyend;
471: lu->MatLUFactorSymbolic = A->ops->lufactorsymbolic;
472: lu->MatDestroy = A->ops->destroy;
473: lu->CleanUpPlapack = PETSC_FALSE;
475: B->spptr = (void*)lu;
476: B->ops->duplicate = MatDuplicate_Plapack;
477: B->ops->view = MatView_Plapack;
478: B->ops->assemblyend = MatAssemblyEnd_Plapack;
479: B->ops->lufactorsymbolic = MatLUFactorSymbolic_Plapack;
480: B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_Plapack;
481: B->ops->destroy = MatDestroy_Plapack;
482:
483: MPI_Comm_size(A->comm,&size);
484: if (size == 1) {
485: PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqdense_plapack_C",
486: "MatConvert_Dense_Plapack",MatConvert_Dense_Plapack);
487: PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_plapack_seqdense_C",
488: "MatConvert_Plapack_Dense",MatConvert_Plapack_Dense);
489: } else {
490: PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_mpidense_plapack_C",
491: "MatConvert_Dense_Plapack",MatConvert_Dense_Plapack);
492: PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_plapack_mpidense_C",
493: "MatConvert_Plapack_Dense",MatConvert_Plapack_Dense);
494: }
495: PetscInfo(0,"Using Plapack for dense LU factorization and solves.\n");
496: PetscObjectChangeTypeName((PetscObject)B,MATPLAPACK);
497: *newmat = B;
498: return(0);
499: }
504: PetscErrorCode MatDuplicate_Plapack(Mat A, MatDuplicateOption op, Mat *M)
505: {
507: Mat_Plapack *lu=(Mat_Plapack *)A->spptr;
510: (*lu->MatDuplicate)(A,op,M);
511: PetscMemcpy((*M)->spptr,lu,sizeof(Mat_Plapack));
512: return(0);
513: }
515: /*MC
516: MATPLAPACK - MATPLAPACK = "plapack" - A matrix type providing direct solvers (LU, Cholesky, and QR)
517: for parallel dense matrices via the external package PLAPACK.
519: If PLAPACK is installed (see the manual for
520: instructions on how to declare the existence of external packages),
521: a matrix type can be constructed which invokes PLAPACK solvers.
522: After calling MatCreate(...,A), simply call MatSetType(A,MATPLAPACK).
524: This matrix inherits from MATSEQDENSE when constructed with a single process communicator,
525: and from MATMPIDENSE otherwise. One can also call MatConvert for an inplace
526: conversion to or from the MATSEQDENSE or MATMPIDENSE type (depending on the communicator size)
527: without data copy.
529: Options Database Keys:
530: + -mat_type plapack - sets the matrix type to "plapack" during a call to MatSetFromOptions()
531: . -mat_plapack_nprows <n> - number of rows in processor partition
532: . -mat_plapack_npcols <n> - number of columns in processor partition
533: . -mat_plapack_nb <n> - block size of template vector
534: . -mat_plapack_nb_alg <n> - algorithmic block size
535: - -mat_plapack_ckerror <n> - error checking flag
537: Level: beginner
539: .seealso: MATDENSE, PCLU, PCCHOLESKY
540: M*/
545: PetscErrorCode MatCreate_Plapack(Mat A)
546: {
548: PetscMPIInt size;
551: MPI_Comm_size(A->comm,&size);
552: if (size == 1) {
553: MatSetType(A,MATSEQDENSE);
554: } else {
555: MatSetType(A,MATMPIDENSE);
556: }
557: MatConvert_Dense_Plapack(A,MATPLAPACK,MAT_REUSE_MATRIX,&A);
558: return(0);
559: }