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

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

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

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

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

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

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

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

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

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




163: /*TEST

165:    test:
166:       nsize: 2
167:       args: -m 3 -n 2
168:       requires: elemental

170:    test:
171:       suffix: 2
172:       nsize: 6
173:       args: -m 2 -n 3
174:       requires: elemental

176:    test:
177:       suffix: 3
178:       nsize: 1
179:       args: -m 2 -n 3
180:       requires: elemental
181:       output_file: output/ex39_1.out

183: TEST*/