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;
 13:   PetscScalar    *v;
 14:   PetscMPIInt    rank,color;
 15:   PetscReal      Cnorm;
 16:   PetscBool      flg,mats_view=PETSC_FALSE;
 17:   MPI_Comm       subcomm;

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

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

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

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

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

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

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

136:   /* Test MatMatTransposeMult(): B = C*C^T */
137:   MatMatTransposeMult(C,C,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&B);
138:   MatScale(C,2.0);
139:   MatMatTransposeMult(C,C,MAT_REUSE_MATRIX,PETSC_DEFAULT,&B);
140:   MatMatTransposeMultEqual(C,C,B,10,&flg);

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

148:   /* Test MatMult() */
149:   MatComputeOperator(C,MATAIJ,&Caij);
150:   MatMultEqual(C,Caij,5,&flg);
152:   MatMultTransposeEqual(C,Caij,5,&flg);

155:   /* Test MatMultAdd() and MatMultTransposeAddEqual() */
156:   MatMultAddEqual(C,Caij,5,&flg);
158:   MatMultTransposeAddEqual(C,Caij,5,&flg);

161:   /* Test MatMatMult() */
162:   PetscOptionsHasName(NULL,NULL,"-test_matmatmult",&flg);
163:   if (flg) {
164:     Mat CC,CCaij;
165:     MatMatMult(C,C,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&CC);
166:     MatMatMult(Caij,Caij,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&CCaij);
167:     MatMultEqual(CC,CCaij,5,&flg);
169:     MatDestroy(&CCaij);
170:     MatDestroy(&CC);
171:   }

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

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

201: /*TEST

203:    build:
204:       requires: scalapack

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

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

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

225: TEST*/