Actual source code: ex39.c


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

  4: #include <petscmat.h>

  6: int main(int argc,char **args)
  7: {
  8:   Mat            Cdense,B,C,Ct;
  9:   Vec            d;
 10:   PetscInt       i,j,m = 5,n,nrows,ncols;
 11:   const PetscInt *rows,*cols;
 12:   IS             isrows,iscols;
 13:   PetscScalar    *v;
 14:   PetscMPIInt    rank,size;
 15:   PetscReal      Cnorm;
 16:   PetscBool      flg,mats_view=PETSC_FALSE;

 18:   PetscInitialize(&argc,&args,(char*)0,help);
 19:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
 20:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 21:   PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL);
 22:   n    = m;
 23:   PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
 24:   PetscOptionsHasName(NULL,NULL,"-mats_view",&mats_view);

 26:   MatCreate(PETSC_COMM_WORLD,&C);
 27:   MatSetSizes(C,m,n,PETSC_DECIDE,PETSC_DECIDE);
 28:   MatSetType(C,MATELEMENTAL);
 29:   MatSetFromOptions(C);
 30:   MatSetUp(C);
 31:   MatGetOwnershipIS(C,&isrows,&iscols);
 32:   ISGetLocalSize(isrows,&nrows);
 33:   ISGetIndices(isrows,&rows);
 34:   ISGetLocalSize(iscols,&ncols);
 35:   ISGetIndices(iscols,&cols);
 36:   PetscMalloc1(nrows*ncols,&v);
 37: #if defined(PETSC_USE_COMPLEX)
 38:   PetscRandom rand;
 39:   PetscScalar rval;
 40:   PetscRandomCreate(PETSC_COMM_WORLD,&rand);
 41:   PetscRandomSetFromOptions(rand);
 42:   for (i=0; i<nrows; i++) {
 43:     for (j=0; j<ncols; j++) {
 44:       PetscRandomGetValue(rand,&rval);
 45:       v[i*ncols+j] = rval;
 46:     }
 47:   }
 48:   PetscRandomDestroy(&rand);
 49: #else
 50:   for (i=0; i<nrows; i++) {
 51:     for (j=0; j<ncols; j++) {
 52:       v[i*ncols+j] = (PetscReal)(10000*rank+100*rows[i]+cols[j]);
 53:     }
 54:   }
 55: #endif
 56:   MatSetValues(C,nrows,rows,ncols,cols,v,INSERT_VALUES);
 57:   ISRestoreIndices(isrows,&rows);
 58:   ISRestoreIndices(iscols,&cols);
 59:   MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
 60:   MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
 61:   ISDestroy(&isrows);
 62:   ISDestroy(&iscols);

 64:   /* Test MatView(), MatDuplicate() and out-of-place MatConvert() */
 65:   MatDuplicate(C,MAT_COPY_VALUES,&B);
 66:   if (mats_view) {
 67:     PetscPrintf(PETSC_COMM_WORLD,"Duplicated C:\n");
 68:     MatView(B,PETSC_VIEWER_STDOUT_WORLD);
 69:   }
 70:   MatDestroy(&B);
 71:   MatConvert(C,MATMPIDENSE,MAT_INITIAL_MATRIX,&Cdense);
 72:   MatMultEqual(C,Cdense,5,&flg);

 75:   /* Test MatNorm() */
 76:   MatNorm(C,NORM_1,&Cnorm);

 78:   /* Test MatTranspose(), MatZeroEntries() and MatGetDiagonal() */
 79:   MatTranspose(C,MAT_INITIAL_MATRIX,&Ct);
 80:   MatConjugate(Ct);
 81:   if (mats_view) {
 82:     PetscPrintf(PETSC_COMM_WORLD,"C's Transpose Conjugate:\n");
 83:     MatView(Ct,PETSC_VIEWER_STDOUT_WORLD);
 84:   }
 85:   MatZeroEntries(Ct);
 86:   VecCreate(PETSC_COMM_WORLD,&d);
 87:   VecSetSizes(d,m>n ? n : m,PETSC_DECIDE);
 88:   VecSetFromOptions(d);
 89:   MatGetDiagonal(C,d);
 90:   if (mats_view) {
 91:     PetscPrintf(PETSC_COMM_WORLD,"Diagonal of C:\n");
 92:     VecView(d,PETSC_VIEWER_STDOUT_WORLD);
 93:   }
 94:   if (m>n) {
 95:     MatDiagonalScale(C,NULL,d);
 96:   } else {
 97:     MatDiagonalScale(C,d,NULL);
 98:   }
 99:   if (mats_view) {
100:     PetscPrintf(PETSC_COMM_WORLD,"Diagonal Scaled C:\n");
101:     MatView(C,PETSC_VIEWER_STDOUT_WORLD);
102:   }

104:   /* Test MatAXPY(), MatAYPX() and in-place MatConvert() */
105:   MatCreate(PETSC_COMM_WORLD,&B);
106:   MatSetSizes(B,m,n,PETSC_DECIDE,PETSC_DECIDE);
107:   MatSetType(B,MATELEMENTAL);
108:   MatSetFromOptions(B);
109:   MatSetUp(B);
110:   MatGetOwnershipIS(B,&isrows,&iscols);
111:   ISGetLocalSize(isrows,&nrows);
112:   ISGetIndices(isrows,&rows);
113:   ISGetLocalSize(iscols,&ncols);
114:   ISGetIndices(iscols,&cols);
115:   for (i=0; i<nrows; i++) {
116:     for (j=0; j<ncols; j++) {
117:       v[i*ncols+j] = (PetscReal)(1000*rows[i]+cols[j]);
118:     }
119:   }
120:   MatSetValues(B,nrows,rows,ncols,cols,v,INSERT_VALUES);
121:   PetscFree(v);
122:   ISRestoreIndices(isrows,&rows);
123:   ISRestoreIndices(iscols,&cols);
124:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
125:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
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);

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

149:   MatDestroy(&Cdense);
150:   PetscFree(v);
151:   MatDestroy(&B);
152:   MatDestroy(&C);
153:   MatDestroy(&Ct);
154:   VecDestroy(&d);
155:   PetscFinalize();
156:   return 0;
157: }

159: /*TEST

161:    test:
162:       nsize: 2
163:       args: -m 3 -n 2
164:       requires: elemental

166:    test:
167:       suffix: 2
168:       nsize: 6
169:       args: -m 2 -n 3
170:       requires: elemental

172:    test:
173:       suffix: 3
174:       nsize: 1
175:       args: -m 2 -n 3
176:       requires: elemental
177:       output_file: output/ex39_1.out

179: TEST*/