Actual source code: matelem.cxx
petsc-3.5.4 2015-05-23
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: elem::Print( *a->emat, "Elemental matrix (cyclic ordering)" );
80: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
81: if (A->factortype == MAT_FACTOR_NONE){
82: Mat Adense;
83: PetscPrintf(PetscObjectComm((PetscObject)viewer),"Elemental matrix (explicit ordering)\n");
84: MatConvert(A,MATDENSE,MAT_INITIAL_MATRIX,&Adense);
85: MatView(Adense,viewer);
86: MatDestroy(&Adense);
87: }
88: } else SETERRQ(PetscObjectComm((PetscObject)viewer),PETSC_ERR_SUP,"Format");
89: } else {
90: /* convert to dense format and call MatView() */
91: Mat Adense;
92: PetscPrintf(PetscObjectComm((PetscObject)viewer),"Elemental matrix (explicit ordering)\n");
93: MatConvert(A,MATDENSE,MAT_INITIAL_MATRIX,&Adense);
94: MatView(Adense,viewer);
95: MatDestroy(&Adense);
96: }
97: return(0);
98: }
102: static PetscErrorCode MatGetInfo_Elemental(Mat A,MatInfoType flag,MatInfo *info)
103: {
104: Mat_Elemental *a = (Mat_Elemental*)A->data;
105: PetscMPIInt rank;
108: MPI_Comm_rank(PetscObjectComm((PetscObject)A),&rank);
110: /* if (!rank) printf(" .........MatGetInfo_Elemental ...\n"); */
111: info->block_size = 1.0;
113: if (flag == MAT_LOCAL) {
114: info->nz_allocated = (double)(*a->emat).AllocatedMemory(); /* locally allocated */
115: info->nz_used = info->nz_allocated;
116: } else if (flag == MAT_GLOBAL_MAX) {
117: //MPI_Allreduce(isend,irecv,5,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)matin));
118: /* see MatGetInfo_MPIAIJ() for getting global info->nz_allocated! */
119: //SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP," MAT_GLOBAL_MAX not written yet");
120: } else if (flag == MAT_GLOBAL_SUM) {
121: //SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP," MAT_GLOBAL_SUM not written yet");
122: info->nz_allocated = (double)(*a->emat).AllocatedMemory(); /* locally allocated */
123: info->nz_used = info->nz_allocated; /* assume Elemental does accurate allocation */
124: //MPI_Allreduce(isend,irecv,1,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)A));
125: //PetscPrintf(PETSC_COMM_SELF," ... [%d] locally allocated %g\n",rank,info->nz_allocated);
126: }
128: info->nz_unneeded = 0.0;
129: info->assemblies = (double)A->num_ass;
130: info->mallocs = 0;
131: info->memory = ((PetscObject)A)->mem;
132: info->fill_ratio_given = 0; /* determined by Elemental */
133: info->fill_ratio_needed = 0;
134: info->factor_mallocs = 0;
135: return(0);
136: }
140: static PetscErrorCode MatSetValues_Elemental(Mat A,PetscInt nr,const PetscInt *rows,PetscInt nc,const PetscInt *cols,const PetscScalar *vals,InsertMode imode)
141: {
143: Mat_Elemental *a = (Mat_Elemental*)A->data;
144: PetscMPIInt rank;
145: PetscInt i,j,rrank,ridx,crank,cidx;
148: MPI_Comm_rank(PetscObjectComm((PetscObject)A),&rank);
150: const elem::Grid &grid = a->emat->Grid();
151: for (i=0; i<nr; i++) {
152: PetscInt erow,ecol,elrow,elcol;
153: if (rows[i] < 0) continue;
154: P2RO(A,0,rows[i],&rrank,&ridx);
155: RO2E(A,0,rrank,ridx,&erow);
156: if (rrank < 0 || ridx < 0 || erow < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Incorrect row translation");
157: for (j=0; j<nc; j++) {
158: if (cols[j] < 0) continue;
159: P2RO(A,1,cols[j],&crank,&cidx);
160: RO2E(A,1,crank,cidx,&ecol);
161: if (crank < 0 || cidx < 0 || ecol < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Incorrect col translation");
162: if (erow % grid.MCSize() != grid.MCRank() || ecol % grid.MRSize() != grid.MRRank()){ /* off-proc entry */
163: if (imode != ADD_VALUES) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only ADD_VALUES to off-processor entry is supported");
164: /* 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); */
165: a->esubmat->Set(0,0, (PetscElemScalar)vals[i*nc+j]);
166: a->interface->Axpy(1.0,*(a->esubmat),erow,ecol);
167: continue;
168: }
169: elrow = erow / grid.MCSize();
170: elcol = ecol / grid.MRSize();
171: switch (imode) {
172: case INSERT_VALUES: a->emat->SetLocal(elrow,elcol,(PetscElemScalar)vals[i*nc+j]); break;
173: case ADD_VALUES: a->emat->UpdateLocal(elrow,elcol,(PetscElemScalar)vals[i*nc+j]); break;
174: default: SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"No support for InsertMode %d",(int)imode);
175: }
176: }
177: }
178: return(0);
179: }
183: static PetscErrorCode MatMult_Elemental(Mat A,Vec X,Vec Y)
184: {
185: Mat_Elemental *a = (Mat_Elemental*)A->data;
186: PetscErrorCode ierr;
187: const PetscElemScalar *x;
188: PetscElemScalar *y;
189: PetscElemScalar one = 1,zero = 0;
192: VecGetArrayRead(X,(const PetscScalar **)&x);
193: VecGetArray(Y,(PetscScalar **)&y);
194: { /* Scoping so that constructor is called before pointer is returned */
195: elem::DistMatrix<PetscElemScalar,elem::VC,elem::STAR> xe, ye;
196: xe.LockedAttach(A->cmap->N,1,*a->grid,0,0,x,A->cmap->n);
197: ye.Attach(A->rmap->N,1,*a->grid,0,0,y,A->rmap->n);
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, ye;
220: xe.LockedAttach(A->rmap->N,1,*a->grid,0,0,x,A->rmap->n);
221: ye.Attach(A->cmap->N,1,*a->grid,0,0,y,A->cmap->n);
222: elem::Gemv(elem::TRANSPOSE,one,*a->emat,xe,zero,ye);
223: }
224: VecRestoreArrayRead(X,(const PetscScalar **)&x);
225: VecRestoreArray(Y,(PetscScalar **)&y);
226: return(0);
227: }
231: static PetscErrorCode MatMultAdd_Elemental(Mat A,Vec X,Vec Y,Vec Z)
232: {
233: Mat_Elemental *a = (Mat_Elemental*)A->data;
234: PetscErrorCode ierr;
235: const PetscElemScalar *x;
236: PetscElemScalar *z;
237: PetscElemScalar one = 1;
240: if (Y != Z) {VecCopy(Y,Z);}
241: VecGetArrayRead(X,(const PetscScalar **)&x);
242: VecGetArray(Z,(PetscScalar **)&z);
243: { /* Scoping so that constructor is called before pointer is returned */
244: elem::DistMatrix<PetscElemScalar,elem::VC,elem::STAR> xe, ze;
245: xe.LockedAttach(A->cmap->N,1,*a->grid,0,0,x,A->cmap->n);
246: ze.Attach(A->rmap->N,1,*a->grid,0,0,z,A->rmap->n);
247: elem::Gemv(elem::NORMAL,one,*a->emat,xe,one,ze);
248: }
249: VecRestoreArrayRead(X,(const PetscScalar **)&x);
250: VecRestoreArray(Z,(PetscScalar **)&z);
251: return(0);
252: }
256: static PetscErrorCode MatMultTransposeAdd_Elemental(Mat A,Vec X,Vec Y,Vec Z)
257: {
258: Mat_Elemental *a = (Mat_Elemental*)A->data;
259: PetscErrorCode ierr;
260: const PetscElemScalar *x;
261: PetscElemScalar *z;
262: PetscElemScalar one = 1;
265: if (Y != Z) {VecCopy(Y,Z);}
266: VecGetArrayRead(X,(const PetscScalar **)&x);
267: VecGetArray(Z,(PetscScalar **)&z);
268: { /* Scoping so that constructor is called before pointer is returned */
269: elem::DistMatrix<PetscElemScalar,elem::VC,elem::STAR> xe, ze;
270: xe.LockedAttach(A->rmap->N,1,*a->grid,0,0,x,A->rmap->n);
271: ze.Attach(A->cmap->N,1,*a->grid,0,0,z,A->cmap->n);
272: elem::Gemv(elem::TRANSPOSE,one,*a->emat,xe,one,ze);
273: }
274: VecRestoreArrayRead(X,(const PetscScalar **)&x);
275: VecRestoreArray(Z,(PetscScalar **)&z);
276: return(0);
277: }
281: static PetscErrorCode MatMatMultNumeric_Elemental(Mat A,Mat B,Mat C)
282: {
283: Mat_Elemental *a = (Mat_Elemental*)A->data;
284: Mat_Elemental *b = (Mat_Elemental*)B->data;
285: Mat_Elemental *c = (Mat_Elemental*)C->data;
286: PetscElemScalar one = 1,zero = 0;
289: { /* Scoping so that constructor is called before pointer is returned */
290: elem::Gemm(elem::NORMAL,elem::NORMAL,one,*a->emat,*b->emat,zero,*c->emat);
291: }
292: C->assembled = PETSC_TRUE;
293: return(0);
294: }
298: static PetscErrorCode MatMatMultSymbolic_Elemental(Mat A,Mat B,PetscReal fill,Mat *C)
299: {
301: Mat Ce;
302: MPI_Comm comm;
305: PetscObjectGetComm((PetscObject)A,&comm);
306: MatCreate(comm,&Ce);
307: MatSetSizes(Ce,A->rmap->n,B->cmap->n,PETSC_DECIDE,PETSC_DECIDE);
308: MatSetType(Ce,MATELEMENTAL);
309: MatSetUp(Ce);
310: *C = Ce;
311: return(0);
312: }
316: static PetscErrorCode MatMatMult_Elemental(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C)
317: {
321: if (scall == MAT_INITIAL_MATRIX){
322: PetscLogEventBegin(MAT_MatMultSymbolic,A,B,0,0);
323: MatMatMultSymbolic_Elemental(A,B,1.0,C);
324: PetscLogEventEnd(MAT_MatMultSymbolic,A,B,0,0);
325: }
326: PetscLogEventBegin(MAT_MatMultNumeric,A,B,0,0);
327: MatMatMultNumeric_Elemental(A,B,*C);
328: PetscLogEventEnd(MAT_MatMultNumeric,A,B,0,0);
329: return(0);
330: }
334: static PetscErrorCode MatMatTransposeMultNumeric_Elemental(Mat A,Mat B,Mat C)
335: {
336: Mat_Elemental *a = (Mat_Elemental*)A->data;
337: Mat_Elemental *b = (Mat_Elemental*)B->data;
338: Mat_Elemental *c = (Mat_Elemental*)C->data;
339: PetscElemScalar one = 1,zero = 0;
342: { /* Scoping so that constructor is called before pointer is returned */
343: elem::Gemm(elem::NORMAL,elem::TRANSPOSE,one,*a->emat,*b->emat,zero,*c->emat);
344: }
345: C->assembled = PETSC_TRUE;
346: return(0);
347: }
351: static PetscErrorCode MatMatTransposeMultSymbolic_Elemental(Mat A,Mat B,PetscReal fill,Mat *C)
352: {
354: Mat Ce;
355: MPI_Comm comm;
358: PetscObjectGetComm((PetscObject)A,&comm);
359: MatCreate(comm,&Ce);
360: MatSetSizes(Ce,A->rmap->n,B->rmap->n,PETSC_DECIDE,PETSC_DECIDE);
361: MatSetType(Ce,MATELEMENTAL);
362: MatSetUp(Ce);
363: *C = Ce;
364: return(0);
365: }
369: static PetscErrorCode MatMatTransposeMult_Elemental(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C)
370: {
374: if (scall == MAT_INITIAL_MATRIX){
375: PetscLogEventBegin(MAT_MatTransposeMultSymbolic,A,B,0,0);
376: MatMatMultSymbolic_Elemental(A,B,1.0,C);
377: PetscLogEventEnd(MAT_MatTransposeMultSymbolic,A,B,0,0);
378: }
379: PetscLogEventBegin(MAT_MatTransposeMultNumeric,A,B,0,0);
380: MatMatTransposeMultNumeric_Elemental(A,B,*C);
381: PetscLogEventEnd(MAT_MatTransposeMultNumeric,A,B,0,0);
382: return(0);
383: }
387: static PetscErrorCode MatGetDiagonal_Elemental(Mat A,Vec D)
388: {
389: PetscInt i,nrows,ncols,nD,rrank,ridx,crank,cidx;
390: Mat_Elemental *a = (Mat_Elemental*)A->data;
391: PetscErrorCode ierr;
392: PetscElemScalar v;
393: MPI_Comm comm;
396: PetscObjectGetComm((PetscObject)A,&comm);
397: MatGetSize(A,&nrows,&ncols);
398: nD = nrows>ncols ? ncols : nrows;
399: for (i=0; i<nD; i++) {
400: PetscInt erow,ecol;
401: P2RO(A,0,i,&rrank,&ridx);
402: RO2E(A,0,rrank,ridx,&erow);
403: if (rrank < 0 || ridx < 0 || erow < 0) SETERRQ(comm,PETSC_ERR_PLIB,"Incorrect row translation");
404: P2RO(A,1,i,&crank,&cidx);
405: RO2E(A,1,crank,cidx,&ecol);
406: if (crank < 0 || cidx < 0 || ecol < 0) SETERRQ(comm,PETSC_ERR_PLIB,"Incorrect col translation");
407: v = a->emat->Get(erow,ecol);
408: VecSetValues(D,1,&i,(PetscScalar*)&v,INSERT_VALUES);
409: }
410: VecAssemblyBegin(D);
411: VecAssemblyEnd(D);
412: return(0);
413: }
417: static PetscErrorCode MatDiagonalScale_Elemental(Mat X,Vec L,Vec R)
418: {
419: Mat_Elemental *x = (Mat_Elemental*)X->data;
420: const PetscElemScalar *d;
421: PetscErrorCode ierr;
424: if (R) {
425: VecGetArrayRead(R,(const PetscScalar **)&d);
426: elem::DistMatrix<PetscElemScalar,elem::VC,elem::STAR> de;
427: de.LockedAttach(X->cmap->N,1,*x->grid,0,0,d,X->cmap->n);
428: elem::DiagonalScale(elem::RIGHT,elem::NORMAL,de,*x->emat);
429: VecRestoreArrayRead(R,(const PetscScalar **)&d);
430: }
431: if (L) {
432: VecGetArrayRead(L,(const PetscScalar **)&d);
433: elem::DistMatrix<PetscElemScalar,elem::VC,elem::STAR> de;
434: de.LockedAttach(X->rmap->N,1,*x->grid,0,0,d,X->rmap->n);
435: elem::DiagonalScale(elem::LEFT,elem::NORMAL,de,*x->emat);
436: VecRestoreArrayRead(L,(const PetscScalar **)&d);
437: }
438: return(0);
439: }
443: static PetscErrorCode MatScale_Elemental(Mat X,PetscScalar a)
444: {
445: Mat_Elemental *x = (Mat_Elemental*)X->data;
448: elem::Scale((PetscElemScalar)a,*x->emat);
449: return(0);
450: }
454: static PetscErrorCode MatAXPY_Elemental(Mat Y,PetscScalar a,Mat X,MatStructure str)
455: {
456: Mat_Elemental *x = (Mat_Elemental*)X->data;
457: Mat_Elemental *y = (Mat_Elemental*)Y->data;
461: elem::Axpy((PetscElemScalar)a,*x->emat,*y->emat);
462: PetscObjectStateIncrease((PetscObject)Y);
463: return(0);
464: }
468: static PetscErrorCode MatCopy_Elemental(Mat A,Mat B,MatStructure str)
469: {
470: Mat_Elemental *a=(Mat_Elemental*)A->data;
471: Mat_Elemental *b=(Mat_Elemental*)B->data;
474: elem::Copy(*a->emat,*b->emat);
475: return(0);
476: }
480: static PetscErrorCode MatDuplicate_Elemental(Mat A,MatDuplicateOption op,Mat *B)
481: {
482: Mat Be;
483: MPI_Comm comm;
484: Mat_Elemental *a=(Mat_Elemental*)A->data;
488: PetscObjectGetComm((PetscObject)A,&comm);
489: MatCreate(comm,&Be);
490: MatSetSizes(Be,A->rmap->n,A->cmap->n,PETSC_DECIDE,PETSC_DECIDE);
491: MatSetType(Be,MATELEMENTAL);
492: MatSetUp(Be);
493: *B = Be;
494: if (op == MAT_COPY_VALUES) {
495: Mat_Elemental *b=(Mat_Elemental*)Be->data;
496: elem::Copy(*a->emat,*b->emat);
497: }
498: Be->assembled = PETSC_TRUE;
499: return(0);
500: }
504: static PetscErrorCode MatTranspose_Elemental(Mat A,MatReuse reuse,Mat *B)
505: {
506: Mat Be = *B;
508: MPI_Comm comm;
509: Mat_Elemental *a = (Mat_Elemental*)A->data, *b;
512: PetscObjectGetComm((PetscObject)A,&comm);
513: /* Only out-of-place supported */
514: if (reuse == MAT_INITIAL_MATRIX){
515: MatCreate(comm,&Be);
516: MatSetSizes(Be,A->cmap->n,A->rmap->n,PETSC_DECIDE,PETSC_DECIDE);
517: MatSetType(Be,MATELEMENTAL);
518: MatSetUp(Be);
519: *B = Be;
520: }
521: b = (Mat_Elemental*)Be->data;
522: elem::Transpose(*a->emat,*b->emat);
523: Be->assembled = PETSC_TRUE;
524: return(0);
525: }
529: static PetscErrorCode MatConjugate_Elemental(Mat A)
530: {
531: Mat_Elemental *a = (Mat_Elemental*)A->data;
534: elem::Conjugate(*a->emat);
535: return(0);
536: }
540: static PetscErrorCode MatHermitianTranspose_Elemental(Mat A,MatReuse reuse,Mat *B)
541: {
542: Mat Be = *B;
544: MPI_Comm comm;
545: Mat_Elemental *a = (Mat_Elemental*)A->data, *b;
548: PetscObjectGetComm((PetscObject)A,&comm);
549: /* Only out-of-place supported */
550: if (reuse == MAT_INITIAL_MATRIX){
551: MatCreate(comm,&Be);
552: MatSetSizes(Be,A->cmap->n,A->rmap->n,PETSC_DECIDE,PETSC_DECIDE);
553: MatSetType(Be,MATELEMENTAL);
554: MatSetUp(Be);
555: *B = Be;
556: }
557: b = (Mat_Elemental*)Be->data;
558: elem::Adjoint(*a->emat,*b->emat);
559: Be->assembled = PETSC_TRUE;
560: return(0);
561: }
565: static PetscErrorCode MatSolve_Elemental(Mat A,Vec B,Vec X)
566: {
567: Mat_Elemental *a = (Mat_Elemental*)A->data;
568: PetscErrorCode ierr;
569: PetscElemScalar *x;
572: VecCopy(B,X);
573: VecGetArray(X,(PetscScalar **)&x);
574: elem::DistMatrix<PetscElemScalar,elem::VC,elem::STAR> xe;
575: xe.Attach(A->rmap->N,1,*a->grid,0,0,x,A->rmap->n);
576: elem::DistMatrix<PetscElemScalar,elem::MC,elem::MR> xer(xe);
577: switch (A->factortype) {
578: case MAT_FACTOR_LU:
579: if ((*a->pivot).AllocatedMemory()) {
580: elem::lu::SolveAfter(elem::NORMAL,*a->emat,*a->pivot,xer);
581: elem::Copy(xer,xe);
582: } else {
583: elem::lu::SolveAfter(elem::NORMAL,*a->emat,xer);
584: elem::Copy(xer,xe);
585: }
586: break;
587: case MAT_FACTOR_CHOLESKY:
588: elem::cholesky::SolveAfter(elem::UPPER,elem::NORMAL,*a->emat,xer);
589: elem::Copy(xer,xe);
590: break;
591: default:
592: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unfactored Matrix or Unsupported MatFactorType");
593: break;
594: }
595: VecRestoreArray(X,(PetscScalar **)&x);
596: return(0);
597: }
601: static PetscErrorCode MatSolveAdd_Elemental(Mat A,Vec B,Vec Y,Vec X)
602: {
603: PetscErrorCode ierr;
606: MatSolve_Elemental(A,B,X);
607: VecAXPY(X,1,Y);
608: return(0);
609: }
613: static PetscErrorCode MatMatSolve_Elemental(Mat A,Mat B,Mat X)
614: {
615: Mat_Elemental *a=(Mat_Elemental*)A->data;
616: Mat_Elemental *b=(Mat_Elemental*)B->data;
617: Mat_Elemental *x=(Mat_Elemental*)X->data;
620: elem::Copy(*b->emat,*x->emat);
621: switch (A->factortype) {
622: case MAT_FACTOR_LU:
623: if ((*a->pivot).AllocatedMemory()) {
624: elem::lu::SolveAfter(elem::NORMAL,*a->emat,*a->pivot,*x->emat);
625: } else {
626: elem::lu::SolveAfter(elem::NORMAL,*a->emat,*x->emat);
627: }
628: break;
629: case MAT_FACTOR_CHOLESKY:
630: elem::cholesky::SolveAfter(elem::UPPER,elem::NORMAL,*a->emat,*x->emat);
631: break;
632: default:
633: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unfactored Matrix or Unsupported MatFactorType");
634: break;
635: }
636: return(0);
637: }
641: static PetscErrorCode MatLUFactor_Elemental(Mat A,IS row,IS col,const MatFactorInfo *info)
642: {
643: Mat_Elemental *a = (Mat_Elemental*)A->data;
646: if (info->dtcol){
647: elem::LU(*a->emat,*a->pivot);
648: } else {
649: elem::LU(*a->emat);
650: }
651: A->factortype = MAT_FACTOR_LU;
652: A->assembled = PETSC_TRUE;
653: return(0);
654: }
658: static PetscErrorCode MatLUFactorNumeric_Elemental(Mat F,Mat A,const MatFactorInfo *info)
659: {
663: MatCopy(A,F,SAME_NONZERO_PATTERN);
664: MatLUFactor_Elemental(F,0,0,info);
665: return(0);
666: }
670: static PetscErrorCode MatLUFactorSymbolic_Elemental(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info)
671: {
673: /* F is create and allocated by MatGetFactor_elemental_petsc(), skip this routine. */
674: return(0);
675: }
679: static PetscErrorCode MatCholeskyFactor_Elemental(Mat A,IS perm,const MatFactorInfo *info)
680: {
681: Mat_Elemental *a = (Mat_Elemental*)A->data;
682: elem::DistMatrix<PetscElemScalar,elem::MC,elem::STAR> d;
685: elem::Cholesky(elem::UPPER,*a->emat);
686: A->factortype = MAT_FACTOR_CHOLESKY;
687: A->assembled = PETSC_TRUE;
688: return(0);
689: }
693: static PetscErrorCode MatCholeskyFactorNumeric_Elemental(Mat F,Mat A,const MatFactorInfo *info)
694: {
698: MatCopy(A,F,SAME_NONZERO_PATTERN);
699: MatCholeskyFactor_Elemental(F,0,info);
700: return(0);
701: }
705: static PetscErrorCode MatCholeskyFactorSymbolic_Elemental(Mat F,Mat A,IS perm,const MatFactorInfo *info)
706: {
708: /* F is create and allocated by MatGetFactor_elemental_petsc(), skip this routine. */
709: return(0);
710: }
714: PetscErrorCode MatFactorGetSolverPackage_elemental_elemental(Mat A,const MatSolverPackage *type)
715: {
717: *type = MATSOLVERELEMENTAL;
718: return(0);
719: }
723: static PetscErrorCode MatGetFactor_elemental_elemental(Mat A,MatFactorType ftype,Mat *F)
724: {
725: Mat B;
729: /* Create the factorization matrix */
730: MatCreate(PetscObjectComm((PetscObject)A),&B);
731: MatSetSizes(B,A->rmap->n,A->cmap->n,PETSC_DECIDE,PETSC_DECIDE);
732: MatSetType(B,MATELEMENTAL);
733: MatSetUp(B);
734: B->factortype = ftype;
735: PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverPackage_C",MatFactorGetSolverPackage_elemental_elemental);
736: *F = B;
737: return(0);
738: }
742: static PetscErrorCode MatNorm_Elemental(Mat A,NormType type,PetscReal *nrm)
743: {
744: Mat_Elemental *a=(Mat_Elemental*)A->data;
747: switch (type){
748: case NORM_1:
749: *nrm = elem::OneNorm(*a->emat);
750: break;
751: case NORM_FROBENIUS:
752: *nrm = elem::FrobeniusNorm(*a->emat);
753: break;
754: case NORM_INFINITY:
755: *nrm = elem::InfinityNorm(*a->emat);
756: break;
757: default:
758: printf("Error: unsupported norm type!\n");
759: }
760: return(0);
761: }
765: static PetscErrorCode MatZeroEntries_Elemental(Mat A)
766: {
767: Mat_Elemental *a=(Mat_Elemental*)A->data;
770: elem::Zero(*a->emat);
771: return(0);
772: }
776: static PetscErrorCode MatGetOwnershipIS_Elemental(Mat A,IS *rows,IS *cols)
777: {
778: Mat_Elemental *a = (Mat_Elemental*)A->data;
780: PetscInt i,m,shift,stride,*idx;
783: if (rows) {
784: m = a->emat->LocalHeight();
785: shift = a->emat->ColShift();
786: stride = a->emat->ColStride();
787: PetscMalloc1(m,&idx);
788: for (i=0; i<m; i++) {
789: PetscInt rank,offset;
790: E2RO(A,0,shift+i*stride,&rank,&offset);
791: RO2P(A,0,rank,offset,&idx[i]);
792: }
793: ISCreateGeneral(PETSC_COMM_SELF,m,idx,PETSC_OWN_POINTER,rows);
794: }
795: if (cols) {
796: m = a->emat->LocalWidth();
797: shift = a->emat->RowShift();
798: stride = a->emat->RowStride();
799: PetscMalloc1(m,&idx);
800: for (i=0; i<m; i++) {
801: PetscInt rank,offset;
802: E2RO(A,1,shift+i*stride,&rank,&offset);
803: RO2P(A,1,rank,offset,&idx[i]);
804: }
805: ISCreateGeneral(PETSC_COMM_SELF,m,idx,PETSC_OWN_POINTER,cols);
806: }
807: return(0);
808: }
812: static PetscErrorCode MatConvert_Elemental_Dense(Mat A,MatType newtype,MatReuse reuse,Mat *B)
813: {
814: Mat Bmpi;
815: Mat_Elemental *a = (Mat_Elemental*)A->data;
816: MPI_Comm comm;
817: PetscErrorCode ierr;
818: PetscInt rrank,ridx,crank,cidx,nrows,ncols,i,j;
819: PetscElemScalar v;
820: PetscBool s1,s2,s3;
823: PetscObjectGetComm((PetscObject)A,&comm);
824: PetscStrcmp(newtype,MATDENSE,&s1);
825: PetscStrcmp(newtype,MATSEQDENSE,&s2);
826: PetscStrcmp(newtype,MATMPIDENSE,&s3);
827: if (!s1 && !s2 && !s3) SETERRQ(comm,PETSC_ERR_SUP,"Unsupported New MatType: must be MATDENSE, MATSEQDENSE or MATMPIDENSE");
828: MatCreate(comm,&Bmpi);
829: MatSetSizes(Bmpi,A->rmap->n,A->cmap->n,PETSC_DECIDE,PETSC_DECIDE);
830: MatSetType(Bmpi,MATDENSE);
831: MatSetUp(Bmpi);
832: MatGetSize(A,&nrows,&ncols);
833: for (i=0; i<nrows; i++) {
834: PetscInt erow,ecol;
835: P2RO(A,0,i,&rrank,&ridx);
836: RO2E(A,0,rrank,ridx,&erow);
837: if (rrank < 0 || ridx < 0 || erow < 0) SETERRQ(comm,PETSC_ERR_PLIB,"Incorrect row translation");
838: for (j=0; j<ncols; j++) {
839: P2RO(A,1,j,&crank,&cidx);
840: RO2E(A,1,crank,cidx,&ecol);
841: if (crank < 0 || cidx < 0 || ecol < 0) SETERRQ(comm,PETSC_ERR_PLIB,"Incorrect col translation");
842: v = a->emat->Get(erow,ecol);
843: MatSetValues(Bmpi,1,&i,1,&j,(PetscScalar *)&v,INSERT_VALUES);
844: }
845: }
846: MatAssemblyBegin(Bmpi,MAT_FINAL_ASSEMBLY);
847: MatAssemblyEnd(Bmpi,MAT_FINAL_ASSEMBLY);
848: if (reuse == MAT_REUSE_MATRIX) {
849: MatHeaderReplace(A,Bmpi);
850: } else {
851: *B = Bmpi;
852: }
853: return(0);
854: }
858: static PetscErrorCode MatDestroy_Elemental(Mat A)
859: {
860: Mat_Elemental *a = (Mat_Elemental*)A->data;
861: PetscErrorCode ierr;
862: Mat_Elemental_Grid *commgrid;
863: PetscBool flg;
864: MPI_Comm icomm;
867: a->interface->Detach();
868: delete a->interface;
869: delete a->esubmat;
870: delete a->emat;
872: elem::mpi::Comm cxxcomm(PetscObjectComm((PetscObject)A));
873: PetscCommDuplicate(cxxcomm.comm,&icomm,NULL);
874: MPI_Attr_get(icomm,Petsc_Elemental_keyval,(void**)&commgrid,(int*)&flg);
875: if (--commgrid->grid_refct == 0) {
876: delete commgrid->grid;
877: PetscFree(commgrid);
878: }
879: PetscCommDestroy(&icomm);
880: PetscObjectComposeFunction((PetscObject)A,"MatGetOwnershipIS_C",NULL);
881: PetscObjectComposeFunction((PetscObject)A,"MatGetFactor_petsc_C",NULL);
882: PetscObjectComposeFunction((PetscObject)A,"MatFactorGetSolverPackage_C",NULL);
883: PetscFree(A->data);
884: return(0);
885: }
889: PetscErrorCode MatSetUp_Elemental(Mat A)
890: {
891: Mat_Elemental *a = (Mat_Elemental*)A->data;
893: PetscMPIInt rsize,csize;
896: PetscLayoutSetUp(A->rmap);
897: PetscLayoutSetUp(A->cmap);
899: a->emat->Resize(A->rmap->N,A->cmap->N);
900: elem::Zero(*a->emat);
902: MPI_Comm_size(A->rmap->comm,&rsize);
903: MPI_Comm_size(A->cmap->comm,&csize);
904: if (csize != rsize) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_INCOMP,"Cannot use row and column communicators of different sizes");
905: a->commsize = rsize;
906: a->mr[0] = A->rmap->N % rsize; if (!a->mr[0]) a->mr[0] = rsize;
907: a->mr[1] = A->cmap->N % csize; if (!a->mr[1]) a->mr[1] = csize;
908: a->m[0] = A->rmap->N / rsize + (a->mr[0] != rsize);
909: a->m[1] = A->cmap->N / csize + (a->mr[1] != csize);
910: return(0);
911: }
915: PetscErrorCode MatAssemblyBegin_Elemental(Mat A, MatAssemblyType type)
916: {
917: Mat_Elemental *a = (Mat_Elemental*)A->data;
920: a->interface->Detach();
921: a->interface->Attach(elem::LOCAL_TO_GLOBAL,*(a->emat));
922: return(0);
923: }
927: PetscErrorCode MatAssemblyEnd_Elemental(Mat A, MatAssemblyType type)
928: {
930: /* Currently does nothing */
931: return(0);
932: }
934: /* -------------------------------------------------------------------*/
935: static struct _MatOps MatOps_Values = {
936: MatSetValues_Elemental,
937: 0,
938: 0,
939: MatMult_Elemental,
940: /* 4*/ MatMultAdd_Elemental,
941: MatMultTranspose_Elemental,
942: MatMultTransposeAdd_Elemental,
943: MatSolve_Elemental,
944: MatSolveAdd_Elemental,
945: 0, //MatSolveTranspose_Elemental,
946: /*10*/ 0, //MatSolveTransposeAdd_Elemental,
947: MatLUFactor_Elemental,
948: MatCholeskyFactor_Elemental,
949: 0,
950: MatTranspose_Elemental,
951: /*15*/ MatGetInfo_Elemental,
952: 0,
953: MatGetDiagonal_Elemental,
954: MatDiagonalScale_Elemental,
955: MatNorm_Elemental,
956: /*20*/ MatAssemblyBegin_Elemental,
957: MatAssemblyEnd_Elemental,
958: 0, //MatSetOption_Elemental,
959: MatZeroEntries_Elemental,
960: /*24*/ 0,
961: MatLUFactorSymbolic_Elemental,
962: MatLUFactorNumeric_Elemental,
963: MatCholeskyFactorSymbolic_Elemental,
964: MatCholeskyFactorNumeric_Elemental,
965: /*29*/ MatSetUp_Elemental,
966: 0,
967: 0,
968: 0,
969: 0,
970: /*34*/ MatDuplicate_Elemental,
971: 0,
972: 0,
973: 0,
974: 0,
975: /*39*/ MatAXPY_Elemental,
976: 0,
977: 0,
978: 0,
979: MatCopy_Elemental,
980: /*44*/ 0,
981: MatScale_Elemental,
982: 0,
983: 0,
984: 0,
985: /*49*/ 0,
986: 0,
987: 0,
988: 0,
989: 0,
990: /*54*/ 0,
991: 0,
992: 0,
993: 0,
994: 0,
995: /*59*/ 0,
996: MatDestroy_Elemental,
997: MatView_Elemental,
998: 0,
999: 0,
1000: /*64*/ 0,
1001: 0,
1002: 0,
1003: 0,
1004: 0,
1005: /*69*/ 0,
1006: 0,
1007: MatConvert_Elemental_Dense,
1008: 0,
1009: 0,
1010: /*74*/ 0,
1011: 0,
1012: 0,
1013: 0,
1014: 0,
1015: /*79*/ 0,
1016: 0,
1017: 0,
1018: 0,
1019: 0,
1020: /*84*/ 0,
1021: 0,
1022: 0,
1023: 0,
1024: 0,
1025: /*89*/ MatMatMult_Elemental,
1026: MatMatMultSymbolic_Elemental,
1027: MatMatMultNumeric_Elemental,
1028: 0,
1029: 0,
1030: /*94*/ 0,
1031: MatMatTransposeMult_Elemental,
1032: MatMatTransposeMultSymbolic_Elemental,
1033: MatMatTransposeMultNumeric_Elemental,
1034: 0,
1035: /*99*/ 0,
1036: 0,
1037: 0,
1038: MatConjugate_Elemental,
1039: 0,
1040: /*104*/0,
1041: 0,
1042: 0,
1043: 0,
1044: 0,
1045: /*109*/MatMatSolve_Elemental,
1046: 0,
1047: 0,
1048: 0,
1049: 0,
1050: /*114*/0,
1051: 0,
1052: 0,
1053: 0,
1054: 0,
1055: /*119*/0,
1056: MatHermitianTranspose_Elemental,
1057: 0,
1058: 0,
1059: 0,
1060: /*124*/0,
1061: 0,
1062: 0,
1063: 0,
1064: 0,
1065: /*129*/0,
1066: 0,
1067: 0,
1068: 0,
1069: 0,
1070: /*134*/0,
1071: 0,
1072: 0,
1073: 0,
1074: 0
1075: };
1077: /*MC
1078: MATELEMENTAL = "elemental" - A matrix type for dense matrices using the Elemental package
1080: Options Database Keys:
1081: + -mat_type elemental - sets the matrix type to "elemental" during a call to MatSetFromOptions()
1082: - -mat_elemental_grid_height - sets Grid Height for 2D cyclic ordering of internal matrix
1084: Level: beginner
1086: .seealso: MATDENSE
1087: M*/
1091: PETSC_EXTERN PetscErrorCode MatCreate_Elemental(Mat A)
1092: {
1093: Mat_Elemental *a;
1094: PetscErrorCode ierr;
1095: PetscBool flg,flg1;
1096: Mat_Elemental_Grid *commgrid;
1097: MPI_Comm icomm;
1098: PetscInt optv1;
1101: PetscElementalInitializePackage();
1102: PetscMemcpy(A->ops,&MatOps_Values,sizeof(struct _MatOps));
1103: A->insertmode = NOT_SET_VALUES;
1105: PetscNewLog(A,&a);
1106: A->data = (void*)a;
1108: /* Set up the elemental matrix */
1109: elem::mpi::Comm cxxcomm(PetscObjectComm((PetscObject)A));
1111: /* Grid needs to be shared between multiple Mats on the same communicator, implement by attribute caching on the MPI_Comm */
1112: if (Petsc_Elemental_keyval == MPI_KEYVAL_INVALID) {
1113: MPI_Keyval_create(MPI_NULL_COPY_FN,MPI_NULL_DELETE_FN,&Petsc_Elemental_keyval,(void*)0);
1114: }
1115: PetscCommDuplicate(cxxcomm.comm,&icomm,NULL);
1116: MPI_Attr_get(icomm,Petsc_Elemental_keyval,(void**)&commgrid,(int*)&flg);
1117: if (!flg) {
1118: PetscNewLog(A,&commgrid);
1120: PetscOptionsBegin(PetscObjectComm((PetscObject)A),((PetscObject)A)->prefix,"Elemental Options","Mat");
1121: /* displayed default grid sizes (CommSize,1) are set by us arbitrarily until elem::Grid() is called */
1122: PetscOptionsInt("-mat_elemental_grid_height","Grid Height","None",elem::mpi::Size(cxxcomm),&optv1,&flg1);
1123: if (flg1) {
1124: if (elem::mpi::Size(cxxcomm) % optv1 != 0) {
1125: SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_INCOMP,"Grid Height %D must evenly divide CommSize %D",optv1,(PetscInt)elem::mpi::Size(cxxcomm));
1126: }
1127: commgrid->grid = new elem::Grid(cxxcomm,optv1); /* use user-provided grid height */
1128: } else {
1129: commgrid->grid = new elem::Grid(cxxcomm); /* use Elemental default grid sizes */
1130: }
1131: commgrid->grid_refct = 1;
1132: MPI_Attr_put(icomm,Petsc_Elemental_keyval,(void*)commgrid);
1133: PetscOptionsEnd();
1134: } else {
1135: commgrid->grid_refct++;
1136: }
1137: PetscCommDestroy(&icomm);
1138: a->grid = commgrid->grid;
1139: a->emat = new elem::DistMatrix<PetscElemScalar>(*a->grid);
1140: a->esubmat = new elem::Matrix<PetscElemScalar>(1,1);
1141: a->interface = new elem::AxpyInterface<PetscElemScalar>;
1142: a->pivot = new elem::DistMatrix<PetscInt,elem::VC,elem::STAR>;
1144: /* build cache for off array entries formed */
1145: a->interface->Attach(elem::LOCAL_TO_GLOBAL,*(a->emat));
1147: PetscObjectComposeFunction((PetscObject)A,"MatGetOwnershipIS_C",MatGetOwnershipIS_Elemental);
1148: PetscObjectComposeFunction((PetscObject)A,"MatGetFactor_elemental_C",MatGetFactor_elemental_elemental);
1150: PetscObjectChangeTypeName((PetscObject)A,MATELEMENTAL);
1151: return(0);
1152: }