Actual source code: matelem.cxx
petsc-3.4.5 2014-06-29
1: #include <../src/mat/impls/elemental/matelemimpl.h> /*I "petscmat.h" I*/
3: /*
4: The variable Petsc_Elemental_keyval is used to indicate an MPI attribute that
5: is attached to a communicator, in this case the attribute is a Mat_Elemental_Grid
6: */
7: static PetscMPIInt Petsc_Elemental_keyval = MPI_KEYVAL_INVALID;
11: /*@C
12: PetscElementalInitializePackage - Initialize Elemental package
14: Logically Collective
16: Level: developer
18: .seealso: MATELEMENTAL, PetscElementalFinalizePackage()
19: @*/
20: PetscErrorCode PetscElementalInitializePackage(void)
21: {
25: if (elem::Initialized()) return(0);
26: { /* We have already initialized MPI, so this song and dance is just to pass these variables (which won't be used by Elemental) through the interface that needs references */
27: int zero = 0;
28: char **nothing = 0;
29: elem::Initialize(zero,nothing);
30: }
31: PetscRegisterFinalize(PetscElementalFinalizePackage);
32: return(0);
33: }
37: /*@C
38: PetscElementalFinalizePackage - Finalize Elemental package
40: Logically Collective
42: Level: developer
44: .seealso: MATELEMENTAL, PetscElementalInitializePackage()
45: @*/
46: PetscErrorCode PetscElementalFinalizePackage(void)
47: {
50: elem::Finalize();
51: return(0);
52: }
56: static PetscErrorCode MatView_Elemental(Mat A,PetscViewer viewer)
57: {
59: Mat_Elemental *a = (Mat_Elemental*)A->data;
60: PetscBool iascii;
63: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
64: if (iascii) {
65: PetscViewerFormat format;
66: PetscViewerGetFormat(viewer,&format);
67: if (format == PETSC_VIEWER_ASCII_INFO) {
68: /* call elemental viewing function */
69: PetscViewerASCIIPrintf(viewer,"Elemental run parameters:\n");
70: PetscViewerASCIIPrintf(viewer," allocated entries=%d\n",(*a->emat).AllocatedMemory());
71: PetscViewerASCIIPrintf(viewer," grid height=%d, grid width=%d\n",(*a->emat).Grid().Height(),(*a->emat).Grid().Width());
72: if (format == PETSC_VIEWER_ASCII_FACTOR_INFO) {
73: /* call elemental viewing function */
74: PetscPrintf(PetscObjectComm((PetscObject)viewer),"test matview_elemental 2\n");
75: }
77: } else if (format == PETSC_VIEWER_DEFAULT) {
78: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
79: PetscObjectPrintClassNamePrefixType((PetscObject)A,viewer,"Matrix Object");
80: a->emat->Print("Elemental matrix (cyclic ordering)");
81: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
82: if (A->factortype == MAT_FACTOR_NONE){
83: Mat Adense;
84: PetscPrintf(PetscObjectComm((PetscObject)viewer),"Elemental matrix (explicit ordering)\n");
85: MatConvert(A,MATDENSE,MAT_INITIAL_MATRIX,&Adense);
86: MatView(Adense,viewer);
87: MatDestroy(&Adense);
88: }
89: } else SETERRQ(PetscObjectComm((PetscObject)viewer),PETSC_ERR_SUP,"Format");
90: } else {
91: /* convert to dense format and call MatView() */
92: Mat Adense;
93: PetscPrintf(PetscObjectComm((PetscObject)viewer),"Elemental matrix (explicit ordering)\n");
94: MatConvert(A,MATDENSE,MAT_INITIAL_MATRIX,&Adense);
95: MatView(Adense,viewer);
96: MatDestroy(&Adense);
97: }
98: return(0);
99: }
103: static PetscErrorCode MatGetInfo_Elemental(Mat A,MatInfoType flag,MatInfo *info)
104: {
105: Mat_Elemental *a = (Mat_Elemental*)A->data;
106: PetscMPIInt rank;
109: MPI_Comm_rank(PetscObjectComm((PetscObject)A),&rank);
111: /* if (!rank) printf(" .........MatGetInfo_Elemental ...\n"); */
112: info->block_size = 1.0;
114: if (flag == MAT_LOCAL) {
115: info->nz_allocated = (double)(*a->emat).AllocatedMemory(); /* locally allocated */
116: info->nz_used = info->nz_allocated;
117: } else if (flag == MAT_GLOBAL_MAX) {
118: //MPI_Allreduce(isend,irecv,5,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)matin));
119: /* see MatGetInfo_MPIAIJ() for getting global info->nz_allocated! */
120: //SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP," MAT_GLOBAL_MAX not written yet");
121: } else if (flag == MAT_GLOBAL_SUM) {
122: //SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP," MAT_GLOBAL_SUM not written yet");
123: info->nz_allocated = (double)(*a->emat).AllocatedMemory(); /* locally allocated */
124: info->nz_used = info->nz_allocated; /* assume Elemental does accurate allocation */
125: //MPI_Allreduce(isend,irecv,1,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)A));
126: //PetscPrintf(PETSC_COMM_SELF," ... [%d] locally allocated %g\n",rank,info->nz_allocated);
127: }
129: info->nz_unneeded = 0.0;
130: info->assemblies = (double)A->num_ass;
131: info->mallocs = 0;
132: info->memory = ((PetscObject)A)->mem;
133: info->fill_ratio_given = 0; /* determined by Elemental */
134: info->fill_ratio_needed = 0;
135: info->factor_mallocs = 0;
136: return(0);
137: }
141: static PetscErrorCode MatSetValues_Elemental(Mat A,PetscInt nr,const PetscInt *rows,PetscInt nc,const PetscInt *cols,const PetscScalar *vals,InsertMode imode)
142: {
144: Mat_Elemental *a = (Mat_Elemental*)A->data;
145: PetscMPIInt rank;
146: PetscInt i,j,rrank,ridx,crank,cidx;
149: MPI_Comm_rank(PetscObjectComm((PetscObject)A),&rank);
151: const elem::Grid &grid = a->emat->Grid();
152: for (i=0; i<nr; i++) {
153: PetscInt erow,ecol,elrow,elcol;
154: if (rows[i] < 0) continue;
155: P2RO(A,0,rows[i],&rrank,&ridx);
156: RO2E(A,0,rrank,ridx,&erow);
157: if (rrank < 0 || ridx < 0 || erow < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Incorrect row translation");
158: for (j=0; j<nc; j++) {
159: if (cols[j] < 0) continue;
160: P2RO(A,1,cols[j],&crank,&cidx);
161: RO2E(A,1,crank,cidx,&ecol);
162: if (crank < 0 || cidx < 0 || ecol < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Incorrect col translation");
163: if (erow % grid.MCSize() != grid.MCRank() || ecol % grid.MRSize() != grid.MRRank()){ /* off-proc entry */
164: if (imode != ADD_VALUES) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only ADD_VALUES to off-processor entry is supported");
165: /* PetscPrintf(PETSC_COMM_SELF,"[%D] add off-proc entry (%D,%D, %g) (%D %D)\n",rank,rows[i],cols[j],*(vals+i*nc),erow,ecol); */
166: a->esubmat->Set(0,0, (PetscElemScalar)vals[i*nc+j]);
167: a->interface->Axpy(1.0,*(a->esubmat),erow,ecol);
168: continue;
169: }
170: elrow = erow / grid.MCSize();
171: elcol = ecol / grid.MRSize();
172: switch (imode) {
173: case INSERT_VALUES: a->emat->SetLocal(elrow,elcol,(PetscElemScalar)vals[i*nc+j]); break;
174: case ADD_VALUES: a->emat->UpdateLocal(elrow,elcol,(PetscElemScalar)vals[i*nc+j]); break;
175: default: SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"No support for InsertMode %d",(int)imode);
176: }
177: }
178: }
179: return(0);
180: }
184: static PetscErrorCode MatMult_Elemental(Mat A,Vec X,Vec Y)
185: {
186: Mat_Elemental *a = (Mat_Elemental*)A->data;
187: PetscErrorCode ierr;
188: const PetscElemScalar *x;
189: PetscElemScalar *y;
190: PetscElemScalar one = 1,zero = 0;
193: VecGetArrayRead(X,(const PetscScalar **)&x);
194: VecGetArray(Y,(PetscScalar **)&y);
195: { /* Scoping so that constructor is called before pointer is returned */
196: elem::DistMatrix<PetscElemScalar,elem::VC,elem::STAR> xe(A->cmap->N,1,0,x,A->cmap->n,*a->grid);
197: elem::DistMatrix<PetscElemScalar,elem::VC,elem::STAR> ye(A->rmap->N,1,0,y,A->rmap->n,*a->grid);
198: elem::Gemv(elem::NORMAL,one,*a->emat,xe,zero,ye);
199: }
200: VecRestoreArrayRead(X,(const PetscScalar **)&x);
201: VecRestoreArray(Y,(PetscScalar **)&y);
202: return(0);
203: }
207: static PetscErrorCode MatMultTranspose_Elemental(Mat A,Vec X,Vec Y)
208: {
209: Mat_Elemental *a = (Mat_Elemental*)A->data;
210: PetscErrorCode ierr;
211: const PetscElemScalar *x;
212: PetscElemScalar *y;
213: PetscElemScalar one = 1,zero = 0;
216: VecGetArrayRead(X,(const PetscScalar **)&x);
217: VecGetArray(Y,(PetscScalar **)&y);
218: { /* Scoping so that constructor is called before pointer is returned */
219: elem::DistMatrix<PetscElemScalar,elem::VC,elem::STAR> xe(A->rmap->N,1,0,x,A->rmap->n,*a->grid);
220: elem::DistMatrix<PetscElemScalar,elem::VC,elem::STAR> ye(A->cmap->N,1,0,y,A->cmap->n,*a->grid);
221: elem::Gemv(elem::TRANSPOSE,one,*a->emat,xe,zero,ye);
222: }
223: VecRestoreArrayRead(X,(const PetscScalar **)&x);
224: VecRestoreArray(Y,(PetscScalar **)&y);
225: return(0);
226: }
230: static PetscErrorCode MatMultAdd_Elemental(Mat A,Vec X,Vec Y,Vec Z)
231: {
232: Mat_Elemental *a = (Mat_Elemental*)A->data;
233: PetscErrorCode ierr;
234: const PetscElemScalar *x;
235: PetscElemScalar *z;
236: PetscElemScalar one = 1;
239: if (Y != Z) {VecCopy(Y,Z);}
240: VecGetArrayRead(X,(const PetscScalar **)&x);
241: VecGetArray(Z,(PetscScalar **)&z);
242: { /* Scoping so that constructor is called before pointer is returned */
243: elem::DistMatrix<PetscElemScalar,elem::VC,elem::STAR> xe(A->cmap->N,1,0,x,A->cmap->n,*a->grid);
244: elem::DistMatrix<PetscElemScalar,elem::VC,elem::STAR> ze(A->rmap->N,1,0,z,A->rmap->n,*a->grid);
245: elem::Gemv(elem::NORMAL,one,*a->emat,xe,one,ze);
246: }
247: VecRestoreArrayRead(X,(const PetscScalar **)&x);
248: VecRestoreArray(Z,(PetscScalar **)&z);
249: return(0);
250: }
254: static PetscErrorCode MatMultTransposeAdd_Elemental(Mat A,Vec X,Vec Y,Vec Z)
255: {
256: Mat_Elemental *a = (Mat_Elemental*)A->data;
257: PetscErrorCode ierr;
258: const PetscElemScalar *x;
259: PetscElemScalar *z;
260: PetscElemScalar one = 1;
263: if (Y != Z) {VecCopy(Y,Z);}
264: VecGetArrayRead(X,(const PetscScalar **)&x);
265: VecGetArray(Z,(PetscScalar **)&z);
266: { /* Scoping so that constructor is called before pointer is returned */
267: elem::DistMatrix<PetscElemScalar,elem::VC,elem::STAR> xe(A->rmap->N,1,0,x,A->rmap->n,*a->grid);
268: elem::DistMatrix<PetscElemScalar,elem::VC,elem::STAR> ze(A->cmap->N,1,0,z,A->cmap->n,*a->grid);
269: elem::Gemv(elem::TRANSPOSE,one,*a->emat,xe,one,ze);
270: }
271: VecRestoreArrayRead(X,(const PetscScalar **)&x);
272: VecRestoreArray(Z,(PetscScalar **)&z);
273: return(0);
274: }
278: static PetscErrorCode MatMatMultNumeric_Elemental(Mat A,Mat B,Mat C)
279: {
280: Mat_Elemental *a = (Mat_Elemental*)A->data;
281: Mat_Elemental *b = (Mat_Elemental*)B->data;
282: Mat_Elemental *c = (Mat_Elemental*)C->data;
283: PetscElemScalar one = 1,zero = 0;
286: { /* Scoping so that constructor is called before pointer is returned */
287: elem::Gemm(elem::NORMAL,elem::NORMAL,one,*a->emat,*b->emat,zero,*c->emat);
288: }
289: C->assembled = PETSC_TRUE;
290: return(0);
291: }
295: static PetscErrorCode MatMatMultSymbolic_Elemental(Mat A,Mat B,PetscReal fill,Mat *C)
296: {
298: Mat Ce;
299: MPI_Comm comm;
302: PetscObjectGetComm((PetscObject)A,&comm);
303: MatCreate(comm,&Ce);
304: MatSetSizes(Ce,A->rmap->n,B->cmap->n,PETSC_DECIDE,PETSC_DECIDE);
305: MatSetType(Ce,MATELEMENTAL);
306: MatSetUp(Ce);
307: *C = Ce;
308: return(0);
309: }
313: static PetscErrorCode MatMatMult_Elemental(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C)
314: {
318: if (scall == MAT_INITIAL_MATRIX){
319: PetscLogEventBegin(MAT_MatMultSymbolic,A,B,0,0);
320: MatMatMultSymbolic_Elemental(A,B,1.0,C);
321: PetscLogEventEnd(MAT_MatMultSymbolic,A,B,0,0);
322: }
323: PetscLogEventBegin(MAT_MatMultNumeric,A,B,0,0);
324: MatMatMultNumeric_Elemental(A,B,*C);
325: PetscLogEventEnd(MAT_MatMultNumeric,A,B,0,0);
326: return(0);
327: }
331: static PetscErrorCode MatMatTransposeMultNumeric_Elemental(Mat A,Mat B,Mat C)
332: {
333: Mat_Elemental *a = (Mat_Elemental*)A->data;
334: Mat_Elemental *b = (Mat_Elemental*)B->data;
335: Mat_Elemental *c = (Mat_Elemental*)C->data;
336: PetscElemScalar one = 1,zero = 0;
339: { /* Scoping so that constructor is called before pointer is returned */
340: elem::Gemm(elem::NORMAL,elem::TRANSPOSE,one,*a->emat,*b->emat,zero,*c->emat);
341: }
342: C->assembled = PETSC_TRUE;
343: return(0);
344: }
348: static PetscErrorCode MatMatTransposeMultSymbolic_Elemental(Mat A,Mat B,PetscReal fill,Mat *C)
349: {
351: Mat Ce;
352: MPI_Comm comm;
355: PetscObjectGetComm((PetscObject)A,&comm);
356: MatCreate(comm,&Ce);
357: MatSetSizes(Ce,A->rmap->n,B->rmap->n,PETSC_DECIDE,PETSC_DECIDE);
358: MatSetType(Ce,MATELEMENTAL);
359: MatSetUp(Ce);
360: *C = Ce;
361: return(0);
362: }
366: static PetscErrorCode MatMatTransposeMult_Elemental(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C)
367: {
371: if (scall == MAT_INITIAL_MATRIX){
372: PetscLogEventBegin(MAT_MatTransposeMultSymbolic,A,B,0,0);
373: MatMatMultSymbolic_Elemental(A,B,1.0,C);
374: PetscLogEventEnd(MAT_MatTransposeMultSymbolic,A,B,0,0);
375: }
376: PetscLogEventBegin(MAT_MatTransposeMultNumeric,A,B,0,0);
377: MatMatTransposeMultNumeric_Elemental(A,B,*C);
378: PetscLogEventEnd(MAT_MatTransposeMultNumeric,A,B,0,0);
379: return(0);
380: }
384: static PetscErrorCode MatGetDiagonal_Elemental(Mat A,Vec D)
385: {
386: PetscInt i,nrows,ncols,nD,rrank,ridx,crank,cidx;
387: Mat_Elemental *a = (Mat_Elemental*)A->data;
388: PetscErrorCode ierr;
389: PetscElemScalar v;
390: MPI_Comm comm;
393: PetscObjectGetComm((PetscObject)A,&comm);
394: MatGetSize(A,&nrows,&ncols);
395: nD = nrows>ncols ? ncols : nrows;
396: for (i=0; i<nD; i++) {
397: PetscInt erow,ecol;
398: P2RO(A,0,i,&rrank,&ridx);
399: RO2E(A,0,rrank,ridx,&erow);
400: if (rrank < 0 || ridx < 0 || erow < 0) SETERRQ(comm,PETSC_ERR_PLIB,"Incorrect row translation");
401: P2RO(A,1,i,&crank,&cidx);
402: RO2E(A,1,crank,cidx,&ecol);
403: if (crank < 0 || cidx < 0 || ecol < 0) SETERRQ(comm,PETSC_ERR_PLIB,"Incorrect col translation");
404: v = a->emat->Get(erow,ecol);
405: VecSetValues(D,1,&i,(PetscScalar*)&v,INSERT_VALUES);
406: }
407: VecAssemblyBegin(D);
408: VecAssemblyEnd(D);
409: return(0);
410: }
414: static PetscErrorCode MatDiagonalScale_Elemental(Mat X,Vec L,Vec R)
415: {
416: Mat_Elemental *x = (Mat_Elemental*)X->data;
417: const PetscElemScalar *d;
418: PetscErrorCode ierr;
421: if (R) {
422: VecGetArrayRead(R,(const PetscScalar **)&d);
423: elem::DistMatrix<PetscElemScalar,elem::VC,elem::STAR> de(X->cmap->N,1,0,d,X->cmap->n,*x->grid);
424: elem::DiagonalScale(elem::RIGHT,elem::NORMAL,de,*x->emat);
425: VecRestoreArrayRead(R,(const PetscScalar **)&d);
426: }
427: if (L) {
428: VecGetArrayRead(L,(const PetscScalar **)&d);
429: elem::DistMatrix<PetscElemScalar,elem::VC,elem::STAR> de(X->rmap->N,1,0,d,X->rmap->n,*x->grid);
430: elem::DiagonalScale(elem::LEFT,elem::NORMAL,de,*x->emat);
431: VecRestoreArrayRead(L,(const PetscScalar **)&d);
432: }
433: return(0);
434: }
438: static PetscErrorCode MatScale_Elemental(Mat X,PetscScalar a)
439: {
440: Mat_Elemental *x = (Mat_Elemental*)X->data;
443: elem::Scal((PetscElemScalar)a,*x->emat);
444: return(0);
445: }
449: static PetscErrorCode MatAXPY_Elemental(Mat Y,PetscScalar a,Mat X,MatStructure str)
450: {
451: Mat_Elemental *x = (Mat_Elemental*)X->data;
452: Mat_Elemental *y = (Mat_Elemental*)Y->data;
455: elem::Axpy((PetscElemScalar)a,*x->emat,*y->emat);
456: return(0);
457: }
461: static PetscErrorCode MatCopy_Elemental(Mat A,Mat B,MatStructure str)
462: {
463: Mat_Elemental *a=(Mat_Elemental*)A->data;
464: Mat_Elemental *b=(Mat_Elemental*)B->data;
467: elem::Copy(*a->emat,*b->emat);
468: return(0);
469: }
473: static PetscErrorCode MatDuplicate_Elemental(Mat A,MatDuplicateOption op,Mat *B)
474: {
475: Mat Be;
476: MPI_Comm comm;
477: Mat_Elemental *a=(Mat_Elemental*)A->data;
481: PetscObjectGetComm((PetscObject)A,&comm);
482: MatCreate(comm,&Be);
483: MatSetSizes(Be,A->rmap->n,A->cmap->n,PETSC_DECIDE,PETSC_DECIDE);
484: MatSetType(Be,MATELEMENTAL);
485: MatSetUp(Be);
486: *B = Be;
487: if (op == MAT_COPY_VALUES) {
488: Mat_Elemental *b=(Mat_Elemental*)Be->data;
489: elem::Copy(*a->emat,*b->emat);
490: }
491: Be->assembled = PETSC_TRUE;
492: return(0);
493: }
497: static PetscErrorCode MatTranspose_Elemental(Mat A,MatReuse reuse,Mat *B)
498: {
499: Mat Be;
501: MPI_Comm comm;
502: Mat_Elemental *a = (Mat_Elemental*)A->data, *b;
505: PetscObjectGetComm((PetscObject)A,&comm);
506: /* Only out-of-place supported */
507: if (reuse == MAT_INITIAL_MATRIX){
508: MatCreate(comm,&Be);
509: MatSetSizes(Be,A->cmap->n,A->rmap->n,PETSC_DECIDE,PETSC_DECIDE);
510: MatSetType(Be,MATELEMENTAL);
511: MatSetUp(Be);
512: *B = Be;
513: }
514: b = (Mat_Elemental*)Be->data;
515: elem::Transpose(*a->emat,*b->emat);
516: Be->assembled = PETSC_TRUE;
517: return(0);
518: }
522: static PetscErrorCode MatConjugate_Elemental(Mat A)
523: {
524: Mat_Elemental *a = (Mat_Elemental*)A->data;
527: elem::Conjugate(*a->emat);
528: return(0);
529: }
533: static PetscErrorCode MatHermitianTranspose_Elemental(Mat A,MatReuse reuse,Mat *B)
534: {
535: Mat Be;
537: MPI_Comm comm;
538: Mat_Elemental *a = (Mat_Elemental*)A->data, *b;
541: PetscObjectGetComm((PetscObject)A,&comm);
542: /* Only out-of-place supported */
543: if (reuse == MAT_INITIAL_MATRIX){
544: MatCreate(comm,&Be);
545: MatSetSizes(Be,A->cmap->n,A->rmap->n,PETSC_DECIDE,PETSC_DECIDE);
546: MatSetType(Be,MATELEMENTAL);
547: MatSetUp(Be);
548: *B = Be;
549: }
550: b = (Mat_Elemental*)Be->data;
551: elem::Adjoint(*a->emat,*b->emat);
552: Be->assembled = PETSC_TRUE;
553: return(0);
554: }
558: static PetscErrorCode MatSolve_Elemental(Mat A,Vec B,Vec X)
559: {
560: Mat_Elemental *a = (Mat_Elemental*)A->data;
561: PetscErrorCode ierr;
562: PetscElemScalar *x;
565: VecCopy(B,X);
566: VecGetArray(X,(PetscScalar **)&x);
567: elem::DistMatrix<PetscElemScalar,elem::VC,elem::STAR> xe(A->rmap->N,1,0,x,A->rmap->n,*a->grid);
568: elem::DistMatrix<PetscElemScalar,elem::MC,elem::MR> xer = xe;
569: switch (A->factortype) {
570: case MAT_FACTOR_LU:
571: if ((*a->pivot).AllocatedMemory()) {
572: elem::lu::SolveAfter(elem::NORMAL,*a->emat,*a->pivot,xer);
573: elem::Copy(xer,xe);
574: } else {
575: elem::lu::SolveAfter(elem::NORMAL,*a->emat,xer);
576: elem::Copy(xer,xe);
577: }
578: break;
579: case MAT_FACTOR_CHOLESKY:
580: elem::cholesky::SolveAfter(elem::UPPER,elem::NORMAL,*a->emat,xer);
581: elem::Copy(xer,xe);
582: break;
583: default:
584: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unfactored Matrix or Unsupported MatFactorType");
585: break;
586: }
587: VecRestoreArray(X,(PetscScalar **)&x);
588: return(0);
589: }
593: static PetscErrorCode MatSolveAdd_Elemental(Mat A,Vec B,Vec Y,Vec X)
594: {
595: PetscErrorCode ierr;
598: MatSolve_Elemental(A,B,X);
599: VecAXPY(X,1,Y);
600: return(0);
601: }
605: static PetscErrorCode MatMatSolve_Elemental(Mat A,Mat B,Mat X)
606: {
607: Mat_Elemental *a=(Mat_Elemental*)A->data;
608: Mat_Elemental *b=(Mat_Elemental*)B->data;
609: Mat_Elemental *x=(Mat_Elemental*)X->data;
612: elem::Copy(*b->emat,*x->emat);
613: switch (A->factortype) {
614: case MAT_FACTOR_LU:
615: if ((*a->pivot).AllocatedMemory()) {
616: elem::lu::SolveAfter(elem::NORMAL,*a->emat,*a->pivot,*x->emat);
617: } else {
618: elem::lu::SolveAfter(elem::NORMAL,*a->emat,*x->emat);
619: }
620: break;
621: case MAT_FACTOR_CHOLESKY:
622: elem::cholesky::SolveAfter(elem::UPPER,elem::NORMAL,*a->emat,*x->emat);
623: break;
624: default:
625: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unfactored Matrix or Unsupported MatFactorType");
626: break;
627: }
628: return(0);
629: }
633: static PetscErrorCode MatLUFactor_Elemental(Mat A,IS row,IS col,const MatFactorInfo *info)
634: {
635: Mat_Elemental *a = (Mat_Elemental*)A->data;
638: if (info->dtcol){
639: elem::LU(*a->emat,*a->pivot);
640: } else {
641: elem::LU(*a->emat);
642: }
643: A->factortype = MAT_FACTOR_LU;
644: A->assembled = PETSC_TRUE;
645: return(0);
646: }
650: static PetscErrorCode MatLUFactorNumeric_Elemental(Mat F,Mat A,const MatFactorInfo *info)
651: {
655: MatCopy(A,F,SAME_NONZERO_PATTERN);
656: MatLUFactor_Elemental(F,0,0,info);
657: return(0);
658: }
662: static PetscErrorCode MatLUFactorSymbolic_Elemental(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info)
663: {
665: /* F is create and allocated by MatGetFactor_elemental_petsc(), skip this routine. */
666: return(0);
667: }
671: static PetscErrorCode MatCholeskyFactor_Elemental(Mat A,IS perm,const MatFactorInfo *info)
672: {
673: Mat_Elemental *a = (Mat_Elemental*)A->data;
674: elem::DistMatrix<PetscElemScalar,elem::MC,elem::STAR> d;
677: elem::Cholesky(elem::UPPER,*a->emat);
678: A->factortype = MAT_FACTOR_CHOLESKY;
679: A->assembled = PETSC_TRUE;
680: return(0);
681: }
685: static PetscErrorCode MatCholeskyFactorNumeric_Elemental(Mat F,Mat A,const MatFactorInfo *info)
686: {
690: MatCopy(A,F,SAME_NONZERO_PATTERN);
691: MatCholeskyFactor_Elemental(F,0,info);
692: return(0);
693: }
697: static PetscErrorCode MatCholeskyFactorSymbolic_Elemental(Mat F,Mat A,IS perm,const MatFactorInfo *info)
698: {
700: /* F is create and allocated by MatGetFactor_elemental_petsc(), skip this routine. */
701: return(0);
702: }
706: PetscErrorCode MatFactorGetSolverPackage_elemental_elemental(Mat A,const MatSolverPackage *type)
707: {
709: *type = MATSOLVERELEMENTAL;
710: return(0);
711: }
715: static PetscErrorCode MatGetFactor_elemental_elemental(Mat A,MatFactorType ftype,Mat *F)
716: {
717: Mat B;
721: /* Create the factorization matrix */
722: MatCreate(PetscObjectComm((PetscObject)A),&B);
723: MatSetSizes(B,A->rmap->n,A->cmap->n,PETSC_DECIDE,PETSC_DECIDE);
724: MatSetType(B,MATELEMENTAL);
725: MatSetUp(B);
726: B->factortype = ftype;
727: PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverPackage_C",MatFactorGetSolverPackage_elemental_elemental);
728: *F = B;
729: return(0);
730: }
734: static PetscErrorCode MatNorm_Elemental(Mat A,NormType type,PetscReal *nrm)
735: {
736: Mat_Elemental *a=(Mat_Elemental*)A->data;
739: switch (type){
740: case NORM_1:
741: *nrm = elem::OneNorm(*a->emat);
742: break;
743: case NORM_FROBENIUS:
744: *nrm = elem::FrobeniusNorm(*a->emat);
745: break;
746: case NORM_INFINITY:
747: *nrm = elem::InfinityNorm(*a->emat);
748: break;
749: default:
750: printf("Error: unsupported norm type!\n");
751: }
752: return(0);
753: }
757: static PetscErrorCode MatZeroEntries_Elemental(Mat A)
758: {
759: Mat_Elemental *a=(Mat_Elemental*)A->data;
762: elem::Zero(*a->emat);
763: return(0);
764: }
768: static PetscErrorCode MatGetOwnershipIS_Elemental(Mat A,IS *rows,IS *cols)
769: {
770: Mat_Elemental *a = (Mat_Elemental*)A->data;
772: PetscInt i,m,shift,stride,*idx;
775: if (rows) {
776: m = a->emat->LocalHeight();
777: shift = a->emat->ColShift();
778: stride = a->emat->ColStride();
779: PetscMalloc(m*sizeof(PetscInt),&idx);
780: for (i=0; i<m; i++) {
781: PetscInt rank,offset;
782: E2RO(A,0,shift+i*stride,&rank,&offset);
783: RO2P(A,0,rank,offset,&idx[i]);
784: }
785: ISCreateGeneral(PETSC_COMM_SELF,m,idx,PETSC_OWN_POINTER,rows);
786: }
787: if (cols) {
788: m = a->emat->LocalWidth();
789: shift = a->emat->RowShift();
790: stride = a->emat->RowStride();
791: PetscMalloc(m*sizeof(PetscInt),&idx);
792: for (i=0; i<m; i++) {
793: PetscInt rank,offset;
794: E2RO(A,1,shift+i*stride,&rank,&offset);
795: RO2P(A,1,rank,offset,&idx[i]);
796: }
797: ISCreateGeneral(PETSC_COMM_SELF,m,idx,PETSC_OWN_POINTER,cols);
798: }
799: return(0);
800: }
804: static PetscErrorCode MatConvert_Elemental_Dense(Mat A,MatType newtype,MatReuse reuse,Mat *B)
805: {
806: Mat Bmpi;
807: Mat_Elemental *a = (Mat_Elemental*)A->data;
808: MPI_Comm comm;
809: PetscErrorCode ierr;
810: PetscInt rrank,ridx,crank,cidx,nrows,ncols,i,j;
811: PetscElemScalar v;
812: PetscBool s1,s2,s3;
815: PetscObjectGetComm((PetscObject)A,&comm);
816: PetscStrcmp(newtype,MATDENSE,&s1);
817: PetscStrcmp(newtype,MATSEQDENSE,&s2);
818: PetscStrcmp(newtype,MATMPIDENSE,&s3);
819: if (!s1 && !s2 && !s3) SETERRQ(comm,PETSC_ERR_SUP,"Unsupported New MatType: must be MATDENSE, MATSEQDENSE or MATMPIDENSE");
820: MatCreate(comm,&Bmpi);
821: MatSetSizes(Bmpi,A->rmap->n,A->cmap->n,PETSC_DECIDE,PETSC_DECIDE);
822: MatSetType(Bmpi,MATDENSE);
823: MatSetUp(Bmpi);
824: MatGetSize(A,&nrows,&ncols);
825: for (i=0; i<nrows; i++) {
826: PetscInt erow,ecol;
827: P2RO(A,0,i,&rrank,&ridx);
828: RO2E(A,0,rrank,ridx,&erow);
829: if (rrank < 0 || ridx < 0 || erow < 0) SETERRQ(comm,PETSC_ERR_PLIB,"Incorrect row translation");
830: for (j=0; j<ncols; j++) {
831: P2RO(A,1,j,&crank,&cidx);
832: RO2E(A,1,crank,cidx,&ecol);
833: if (crank < 0 || cidx < 0 || ecol < 0) SETERRQ(comm,PETSC_ERR_PLIB,"Incorrect col translation");
834: v = a->emat->Get(erow,ecol);
835: MatSetValues(Bmpi,1,&i,1,&j,(PetscScalar *)&v,INSERT_VALUES);
836: }
837: }
838: MatAssemblyBegin(Bmpi,MAT_FINAL_ASSEMBLY);
839: MatAssemblyEnd(Bmpi,MAT_FINAL_ASSEMBLY);
840: if (reuse == MAT_REUSE_MATRIX) {
841: MatHeaderReplace(A,Bmpi);
842: } else {
843: *B = Bmpi;
844: }
845: return(0);
846: }
850: static PetscErrorCode MatDestroy_Elemental(Mat A)
851: {
852: Mat_Elemental *a = (Mat_Elemental*)A->data;
853: PetscErrorCode ierr;
854: Mat_Elemental_Grid *commgrid;
855: PetscBool flg;
856: MPI_Comm icomm;
859: a->interface->Detach();
860: delete a->interface;
861: delete a->esubmat;
862: delete a->emat;
864: elem::mpi::Comm cxxcomm(PetscObjectComm((PetscObject)A));
865: PetscCommDuplicate(cxxcomm,&icomm,NULL);
866: MPI_Attr_get(icomm,Petsc_Elemental_keyval,(void**)&commgrid,(int*)&flg);
867: if (--commgrid->grid_refct == 0) {
868: delete commgrid->grid;
869: PetscFree(commgrid);
870: }
871: PetscCommDestroy(&icomm);
872: PetscObjectComposeFunction((PetscObject)A,"MatGetOwnershipIS_C",NULL);
873: PetscObjectComposeFunction((PetscObject)A,"MatGetFactor_petsc_C",NULL);
874: PetscObjectComposeFunction((PetscObject)A,"MatFactorGetSolverPackage_C",NULL);
875: PetscFree(A->data);
876: return(0);
877: }
881: PetscErrorCode MatSetUp_Elemental(Mat A)
882: {
883: Mat_Elemental *a = (Mat_Elemental*)A->data;
885: PetscMPIInt rsize,csize;
888: PetscLayoutSetUp(A->rmap);
889: PetscLayoutSetUp(A->cmap);
891: a->emat->ResizeTo(A->rmap->N,A->cmap->N);
892: elem::Zero(*a->emat);
894: MPI_Comm_size(A->rmap->comm,&rsize);
895: MPI_Comm_size(A->cmap->comm,&csize);
896: if (csize != rsize) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_INCOMP,"Cannot use row and column communicators of different sizes");
897: a->commsize = rsize;
898: a->mr[0] = A->rmap->N % rsize; if (!a->mr[0]) a->mr[0] = rsize;
899: a->mr[1] = A->cmap->N % csize; if (!a->mr[1]) a->mr[1] = csize;
900: a->m[0] = A->rmap->N / rsize + (a->mr[0] != rsize);
901: a->m[1] = A->cmap->N / csize + (a->mr[1] != csize);
902: return(0);
903: }
907: PetscErrorCode MatAssemblyBegin_Elemental(Mat A, MatAssemblyType type)
908: {
909: Mat_Elemental *a = (Mat_Elemental*)A->data;
912: a->interface->Detach();
913: a->interface->Attach(elem::LOCAL_TO_GLOBAL,*(a->emat));
914: return(0);
915: }
919: PetscErrorCode MatAssemblyEnd_Elemental(Mat A, MatAssemblyType type)
920: {
922: /* Currently does nothing */
923: return(0);
924: }
926: /* -------------------------------------------------------------------*/
927: static struct _MatOps MatOps_Values = {
928: MatSetValues_Elemental,
929: 0,
930: 0,
931: MatMult_Elemental,
932: /* 4*/ MatMultAdd_Elemental,
933: MatMultTranspose_Elemental,
934: MatMultTransposeAdd_Elemental,
935: MatSolve_Elemental,
936: MatSolveAdd_Elemental,
937: 0, //MatSolveTranspose_Elemental,
938: /*10*/ 0, //MatSolveTransposeAdd_Elemental,
939: MatLUFactor_Elemental,
940: MatCholeskyFactor_Elemental,
941: 0,
942: MatTranspose_Elemental,
943: /*15*/ MatGetInfo_Elemental,
944: 0,
945: MatGetDiagonal_Elemental,
946: MatDiagonalScale_Elemental,
947: MatNorm_Elemental,
948: /*20*/ MatAssemblyBegin_Elemental,
949: MatAssemblyEnd_Elemental,
950: 0, //MatSetOption_Elemental,
951: MatZeroEntries_Elemental,
952: /*24*/ 0,
953: MatLUFactorSymbolic_Elemental,
954: MatLUFactorNumeric_Elemental,
955: MatCholeskyFactorSymbolic_Elemental,
956: MatCholeskyFactorNumeric_Elemental,
957: /*29*/ MatSetUp_Elemental,
958: 0,
959: 0,
960: 0,
961: 0,
962: /*34*/ MatDuplicate_Elemental,
963: 0,
964: 0,
965: 0,
966: 0,
967: /*39*/ MatAXPY_Elemental,
968: 0,
969: 0,
970: 0,
971: MatCopy_Elemental,
972: /*44*/ 0,
973: MatScale_Elemental,
974: 0,
975: 0,
976: 0,
977: /*49*/ 0,
978: 0,
979: 0,
980: 0,
981: 0,
982: /*54*/ 0,
983: 0,
984: 0,
985: 0,
986: 0,
987: /*59*/ 0,
988: MatDestroy_Elemental,
989: MatView_Elemental,
990: 0,
991: 0,
992: /*64*/ 0,
993: 0,
994: 0,
995: 0,
996: 0,
997: /*69*/ 0,
998: 0,
999: MatConvert_Elemental_Dense,
1000: 0,
1001: 0,
1002: /*74*/ 0,
1003: 0,
1004: 0,
1005: 0,
1006: 0,
1007: /*79*/ 0,
1008: 0,
1009: 0,
1010: 0,
1011: 0,
1012: /*84*/ 0,
1013: 0,
1014: 0,
1015: 0,
1016: 0,
1017: /*89*/ MatMatMult_Elemental,
1018: MatMatMultSymbolic_Elemental,
1019: MatMatMultNumeric_Elemental,
1020: 0,
1021: 0,
1022: /*94*/ 0,
1023: MatMatTransposeMult_Elemental,
1024: MatMatTransposeMultSymbolic_Elemental,
1025: MatMatTransposeMultNumeric_Elemental,
1026: 0,
1027: /*99*/ 0,
1028: 0,
1029: 0,
1030: MatConjugate_Elemental,
1031: 0,
1032: /*104*/0,
1033: 0,
1034: 0,
1035: 0,
1036: 0,
1037: /*109*/MatMatSolve_Elemental,
1038: 0,
1039: 0,
1040: 0,
1041: 0,
1042: /*114*/0,
1043: 0,
1044: 0,
1045: 0,
1046: 0,
1047: /*119*/0,
1048: MatHermitianTranspose_Elemental,
1049: 0,
1050: 0,
1051: 0,
1052: /*124*/0,
1053: 0,
1054: 0,
1055: 0,
1056: 0,
1057: /*129*/0,
1058: 0,
1059: 0,
1060: 0,
1061: 0,
1062: /*134*/0,
1063: 0,
1064: 0,
1065: 0,
1066: 0
1067: };
1069: /*MC
1070: MATELEMENTAL = "elemental" - A matrix type for dense matrices using the Elemental package
1072: Options Database Keys:
1073: . -mat_type elemental - sets the matrix type to "elemental" during a call to MatSetFromOptions()
1074: . -mat_elemental_grid_height - sets Grid Height
1075: . -mat_elemental_grid_width - sets Grid Width
1077: Level: beginner
1079: .seealso: MATDENSE
1080: M*/
1084: PETSC_EXTERN PetscErrorCode MatCreate_Elemental(Mat A)
1085: {
1086: Mat_Elemental *a;
1087: PetscErrorCode ierr;
1088: PetscBool flg,flg1,flg2;
1089: Mat_Elemental_Grid *commgrid;
1090: MPI_Comm icomm;
1091: PetscInt optv1,optv2;
1094: PetscElementalInitializePackage();
1095: PetscMemcpy(A->ops,&MatOps_Values,sizeof(struct _MatOps));
1096: A->insertmode = NOT_SET_VALUES;
1098: PetscNewLog(A,Mat_Elemental,&a);
1099: A->data = (void*)a;
1101: /* Set up the elemental matrix */
1102: elem::mpi::Comm cxxcomm(PetscObjectComm((PetscObject)A));
1104: /* Grid needs to be shared between multiple Mats on the same communicator, implement by attribute caching on the MPI_Comm */
1105: if (Petsc_Elemental_keyval == MPI_KEYVAL_INVALID) {
1106: MPI_Keyval_create(MPI_NULL_COPY_FN,MPI_NULL_DELETE_FN,&Petsc_Elemental_keyval,(void*)0);
1107: }
1108: PetscCommDuplicate(cxxcomm,&icomm,NULL);
1109: MPI_Attr_get(icomm,Petsc_Elemental_keyval,(void**)&commgrid,(int*)&flg);
1110: if (!flg) {
1111: PetscNewLog(A,Mat_Elemental_Grid,&commgrid);
1113: PetscOptionsBegin(PetscObjectComm((PetscObject)A),((PetscObject)A)->prefix,"Elemental Options","Mat");
1114: /* displayed default grid sizes (CommSize,1) are set by us arbitrarily until elem::Grid() is called */
1115: PetscOptionsInt("-mat_elemental_grid_height","Grid Height","None",elem::mpi::CommSize(cxxcomm),&optv1,&flg1);
1116: PetscOptionsInt("-mat_elemental_grid_width","Grid Width","None",1,&optv2,&flg2);
1117: if (flg1 || flg2) {
1118: if (optv1*optv2 != elem::mpi::CommSize(cxxcomm)) {
1119: SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_INCOMP,"Grid Height times Grid Width must equal CommSize");
1120: }
1121: commgrid->grid = new elem::Grid(cxxcomm,optv1,optv2); /* use user-provided grid sizes */
1122: } else {
1123: commgrid->grid = new elem::Grid(cxxcomm); /* use Elemental default grid sizes */
1124: }
1125: commgrid->grid_refct = 1;
1126: MPI_Attr_put(icomm,Petsc_Elemental_keyval,(void*)commgrid);
1127: PetscOptionsEnd();
1128: } else {
1129: commgrid->grid_refct++;
1130: }
1131: PetscCommDestroy(&icomm);
1132: a->grid = commgrid->grid;
1133: a->emat = new elem::DistMatrix<PetscElemScalar>(*a->grid);
1134: a->esubmat = new elem::Matrix<PetscElemScalar>(1,1);
1135: a->interface = new elem::AxpyInterface<PetscElemScalar>;
1136: a->pivot = new elem::DistMatrix<PetscInt,elem::VC,elem::STAR>;
1138: /* build cache for off array entries formed */
1139: a->interface->Attach(elem::LOCAL_TO_GLOBAL,*(a->emat));
1141: PetscObjectComposeFunction((PetscObject)A,"MatGetOwnershipIS_C",MatGetOwnershipIS_Elemental);
1142: PetscObjectComposeFunction((PetscObject)A,"MatGetFactor_elemental_C",MatGetFactor_elemental_elemental);
1144: PetscObjectChangeTypeName((PetscObject)A,MATELEMENTAL);
1145: return(0);
1146: }