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