Actual source code: ex242.c


  2: static char help[] = "Tests ScaLAPACK interface.\n\n";

  4: #include <petscmat.h>

  6: int main(int argc,char **args)
  7: {
  8:   Mat            Cdense,Caij,B,C,Ct,Asub;
  9:   Vec            d;
 10:   PetscInt       i,j,M = 5,N,mb = 2,nb,nrows,ncols,mloc,nloc;
 11:   const PetscInt *rows,*cols;
 12:   IS             isrows,iscols;
 14:   PetscScalar    *v;
 15:   PetscMPIInt    rank,color;
 16:   PetscReal      Cnorm;
 17:   PetscBool      flg,mats_view=PETSC_FALSE;
 18:   MPI_Comm       subcomm;

 20:   PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
 21:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
 22:   PetscOptionsGetInt(NULL,NULL,"-M",&M,NULL);
 23:   N    = M;
 24:   PetscOptionsGetInt(NULL,NULL,"-N",&N,NULL);
 25:   PetscOptionsGetInt(NULL,NULL,"-mb",&mb,NULL);
 26:   nb   = mb;
 27:   PetscOptionsGetInt(NULL,NULL,"-nb",&nb,NULL);
 28:   PetscOptionsHasName(NULL,NULL,"-mats_view",&mats_view);

 30:   MatCreate(PETSC_COMM_WORLD,&C);
 31:   MatSetType(C,MATSCALAPACK);
 32:   mloc = PETSC_DECIDE;
 33:   PetscSplitOwnershipEqual(PETSC_COMM_WORLD,&mloc,&M);
 34:   nloc = PETSC_DECIDE;
 35:   PetscSplitOwnershipEqual(PETSC_COMM_WORLD,&nloc,&N);
 36:   MatSetSizes(C,mloc,nloc,M,N);
 37:   MatScaLAPACKSetBlockSizes(C,mb,nb);
 38:   MatSetFromOptions(C);
 39:   MatSetUp(C);
 40:   /*MatCreateScaLAPACK(PETSC_COMM_WORLD,mb,nb,M,N,0,0,&C); */

 42:   MatGetOwnershipIS(C,&isrows,&iscols);
 43:   ISGetLocalSize(isrows,&nrows);
 44:   ISGetIndices(isrows,&rows);
 45:   ISGetLocalSize(iscols,&ncols);
 46:   ISGetIndices(iscols,&cols);
 47:   PetscMalloc1(nrows*ncols,&v);
 48:   for (i=0;i<nrows;i++) {
 49:     for (j=0;j<ncols;j++) v[i*ncols+j] = (PetscReal)(rows[i]+1+(cols[j]+1)*0.01);
 50:   }
 51:   MatSetValues(C,nrows,rows,ncols,cols,v,INSERT_VALUES);
 52:   PetscFree(v);
 53:   ISRestoreIndices(isrows,&rows);
 54:   ISRestoreIndices(iscols,&cols);
 55:   MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
 56:   MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
 57:   ISDestroy(&isrows);
 58:   ISDestroy(&iscols);

 60:   /* Test MatView(), MatDuplicate() and out-of-place MatConvert() */
 61:   MatDuplicate(C,MAT_COPY_VALUES,&B);
 62:   if (mats_view) {
 63:     PetscPrintf(PETSC_COMM_WORLD,"Duplicated C:\n");
 64:     MatView(B,PETSC_VIEWER_STDOUT_WORLD);
 65:   }
 66:   MatDestroy(&B);
 67:   MatConvert(C,MATDENSE,MAT_INITIAL_MATRIX,&Cdense);
 68:   MatMultEqual(C,Cdense,5,&flg);
 69:   if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Check fails: Cdense != C");

 71:   /* Test MatNorm() */
 72:   MatNorm(C,NORM_1,&Cnorm);

 74:   /* Test MatTranspose(), MatZeroEntries() and MatGetDiagonal() */
 75:   MatTranspose(C,MAT_INITIAL_MATRIX,&Ct);
 76:   MatConjugate(Ct);
 77:   if (mats_view) {
 78:     PetscPrintf(PETSC_COMM_WORLD,"C's Transpose Conjugate:\n");
 79:     MatView(Ct,PETSC_VIEWER_STDOUT_WORLD);
 80:   }
 81:   MatZeroEntries(Ct);
 82:   if (M>N) { MatCreateVecs(C,&d,NULL); }
 83:   else { MatCreateVecs(C,NULL,&d); }
 84:   MatGetDiagonal(C,d);
 85:   if (mats_view) {
 86:     PetscPrintf(PETSC_COMM_WORLD,"Diagonal of C:\n");
 87:     VecView(d,PETSC_VIEWER_STDOUT_WORLD);
 88:   }
 89:   if (M>N) {
 90:     MatDiagonalScale(C,NULL,d);
 91:   } else {
 92:     MatDiagonalScale(C,d,NULL);
 93:   }
 94:   if (mats_view) {
 95:     PetscPrintf(PETSC_COMM_WORLD,"Diagonal Scaled C:\n");
 96:     MatView(C,PETSC_VIEWER_STDOUT_WORLD);
 97:   }

 99:   /* Test MatAXPY(), MatAYPX() and in-place MatConvert() */
100:   MatCreate(PETSC_COMM_WORLD,&B);
101:   MatSetType(B,MATSCALAPACK);
102:   MatSetSizes(B,mloc,nloc,PETSC_DECIDE,PETSC_DECIDE);
103:   MatScaLAPACKSetBlockSizes(B,mb,nb);
104:   MatSetFromOptions(B);
105:   MatSetUp(B);
106:   /* MatCreateScaLAPACK(PETSC_COMM_WORLD,mb,nb,M,N,0,0,&B); */
107:   MatGetOwnershipIS(B,&isrows,&iscols);
108:   ISGetLocalSize(isrows,&nrows);
109:   ISGetIndices(isrows,&rows);
110:   ISGetLocalSize(iscols,&ncols);
111:   ISGetIndices(iscols,&cols);
112:   PetscMalloc1(nrows*ncols,&v);
113:   for (i=0;i<nrows;i++) {
114:     for (j=0;j<ncols;j++) v[i*ncols+j] = (PetscReal)(1000*rows[i]+cols[j]);
115:   }
116:   MatSetValues(B,nrows,rows,ncols,cols,v,INSERT_VALUES);
117:   PetscFree(v);
118:   ISRestoreIndices(isrows,&rows);
119:   ISRestoreIndices(iscols,&cols);
120:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
121:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
122:   if (mats_view) {
123:     PetscPrintf(PETSC_COMM_WORLD,"B:\n");
124:     MatView(B,PETSC_VIEWER_STDOUT_WORLD);
125:   }
126:   MatAXPY(B,2.5,C,SAME_NONZERO_PATTERN);
127:   MatAYPX(B,3.75,C,SAME_NONZERO_PATTERN);
128:   MatConvert(B,MATDENSE,MAT_INPLACE_MATRIX,&B);
129:   if (mats_view) {
130:     PetscPrintf(PETSC_COMM_WORLD,"B after MatAXPY and MatAYPX:\n");
131:     MatView(B,PETSC_VIEWER_STDOUT_WORLD);
132:   }
133:   ISDestroy(&isrows);
134:   ISDestroy(&iscols);
135:   MatDestroy(&B);

137:   /* Test MatMatTransposeMult(): B = C*C^T */
138:   MatMatTransposeMult(C,C,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&B);
139:   MatScale(C,2.0);
140:   MatMatTransposeMult(C,C,MAT_REUSE_MATRIX,PETSC_DEFAULT,&B);
141:   MatMatTransposeMultEqual(C,C,B,10,&flg);
142:   if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Check fails: B != C*C^T");

144:   if (mats_view) {
145:     PetscPrintf(PETSC_COMM_WORLD,"C MatMatTransposeMult C:\n");
146:     MatView(B,PETSC_VIEWER_STDOUT_WORLD);
147:   }

149:   /* Test MatMult() */
150:   MatComputeOperator(C,MATAIJ,&Caij);
151:   MatMultEqual(C,Caij,5,&flg);
152:   if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_NOTSAMETYPE,"C != Caij. MatMultEqual() fails");
153:   MatMultTransposeEqual(C,Caij,5,&flg);
154:   if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_NOTSAMETYPE,"C != Caij. MatMultTransposeEqual() fails");

156:   /* Test MatMultAdd() and MatMultTransposeAddEqual() */
157:   MatMultAddEqual(C,Caij,5,&flg);
158:   if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_NOTSAMETYPE,"C != Caij. MatMultAddEqual() fails");
159:   MatMultTransposeAddEqual(C,Caij,5,&flg);
160:   if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_NOTSAMETYPE,"C != Caij. MatMultTransposeAddEqual() fails");

162:   /* Test MatMatMult() */
163:   PetscOptionsHasName(NULL,NULL,"-test_matmatmult",&flg);
164:   if (flg) {
165:     Mat CC,CCaij;
166:     MatMatMult(C,C,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&CC);
167:     MatMatMult(Caij,Caij,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&CCaij);
168:     MatMultEqual(CC,CCaij,5,&flg);
169:     if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_NOTSAMETYPE,"CC != CCaij. MatMatMult() fails");
170:     MatDestroy(&CCaij);
171:     MatDestroy(&CC);
172:   }

174:   /* Test MatCreate() on subcomm */
175:   color = rank%2;
176:   MPI_Comm_split(PETSC_COMM_WORLD,color,0,&subcomm);
177:   if (color==0) {
178:     MatCreate(subcomm,&Asub);
179:     MatSetType(Asub,MATSCALAPACK);
180:     mloc = PETSC_DECIDE;
181:     PetscSplitOwnershipEqual(subcomm,&mloc,&M);
182:     nloc = PETSC_DECIDE;
183:     PetscSplitOwnershipEqual(subcomm,&nloc,&N);
184:     MatSetSizes(Asub,mloc,nloc,M,N);
185:     MatScaLAPACKSetBlockSizes(Asub,mb,nb);
186:     MatSetFromOptions(Asub);
187:     MatSetUp(Asub);
188:     MatDestroy(&Asub);
189:   }

191:   MatDestroy(&Cdense);
192:   MatDestroy(&Caij);
193:   MatDestroy(&B);
194:   MatDestroy(&C);
195:   MatDestroy(&Ct);
196:   VecDestroy(&d);
197:   MPI_Comm_free(&subcomm);
198:   PetscFinalize();
199:   return ierr;
200: }

202: /*TEST

204:    build:
205:       requires: scalapack

207:    test:
208:       nsize: 2
209:       args: -mb 5 -nb 5 -M 12 -N 10
210:       requires: scalapack

212:    test:
213:       suffix: 2
214:       nsize: 6
215:       args: -mb 8 -nb 6 -M 20 -N 50
216:       requires: scalapack
217:       output_file: output/ex242_1.out

219:    test:
220:       suffix: 3
221:       nsize: 3
222:       args: -mb 2 -nb 2 -M 20 -N 20 -test_matmatmult
223:       requires: scalapack
224:       output_file: output/ex242_1.out

226: TEST*/