1 2 static char help[] = "Tests Elemental interface.\n\n"; 3 4 #include <petscmat.h> 5 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; 17 18 PetscFunctionBeginUser; 19 PetscCall(PetscInitialize(&argc,&args,(char*)0,help)); 20 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank)); 21 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 22 PetscCall(PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL)); 23 n = m; 24 PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL)); 25 PetscCall(PetscOptionsHasName(NULL,NULL,"-mats_view",&mats_view)); 26 27 PetscCall(MatCreate(PETSC_COMM_WORLD,&C)); 28 PetscCall(MatSetSizes(C,m,n,PETSC_DECIDE,PETSC_DECIDE)); 29 PetscCall(MatSetType(C,MATELEMENTAL)); 30 PetscCall(MatSetFromOptions(C)); 31 PetscCall(MatSetUp(C)); 32 PetscCall(MatGetOwnershipIS(C,&isrows,&iscols)); 33 PetscCall(ISGetLocalSize(isrows,&nrows)); 34 PetscCall(ISGetIndices(isrows,&rows)); 35 PetscCall(ISGetLocalSize(iscols,&ncols)); 36 PetscCall(ISGetIndices(iscols,&cols)); 37 PetscCall(PetscMalloc1(nrows*ncols,&v)); 38 #if defined(PETSC_USE_COMPLEX) 39 PetscRandom rand; 40 PetscScalar rval; 41 PetscCall(PetscRandomCreate(PETSC_COMM_WORLD,&rand)); 42 PetscCall(PetscRandomSetFromOptions(rand)); 43 for (i=0; i<nrows; i++) { 44 for (j=0; j<ncols; j++) { 45 PetscCall(PetscRandomGetValue(rand,&rval)); 46 v[i*ncols+j] = rval; 47 } 48 } 49 PetscCall(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 PetscCall(MatSetValues(C,nrows,rows,ncols,cols,v,INSERT_VALUES)); 58 PetscCall(ISRestoreIndices(isrows,&rows)); 59 PetscCall(ISRestoreIndices(iscols,&cols)); 60 PetscCall(MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY)); 61 PetscCall(MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY)); 62 PetscCall(ISDestroy(&isrows)); 63 PetscCall(ISDestroy(&iscols)); 64 65 /* Test MatView(), MatDuplicate() and out-of-place MatConvert() */ 66 PetscCall(MatDuplicate(C,MAT_COPY_VALUES,&B)); 67 if (mats_view) { 68 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Duplicated C:\n")); 69 PetscCall(MatView(B,PETSC_VIEWER_STDOUT_WORLD)); 70 } 71 PetscCall(MatDestroy(&B)); 72 PetscCall(MatConvert(C,MATMPIDENSE,MAT_INITIAL_MATRIX,&Cdense)); 73 PetscCall(MatMultEqual(C,Cdense,5,&flg)); 74 PetscCheck(flg,PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMETYPE,"Cdense != C. MatConvert() fails"); 75 76 /* Test MatNorm() */ 77 PetscCall(MatNorm(C,NORM_1,&Cnorm)); 78 79 /* Test MatTranspose(), MatZeroEntries() and MatGetDiagonal() */ 80 PetscCall(MatTranspose(C,MAT_INITIAL_MATRIX,&Ct)); 81 PetscCall(MatConjugate(Ct)); 82 if (mats_view) { 83 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"C's Transpose Conjugate:\n")); 84 PetscCall(MatView(Ct,PETSC_VIEWER_STDOUT_WORLD)); 85 } 86 PetscCall(MatZeroEntries(Ct)); 87 PetscCall(VecCreate(PETSC_COMM_WORLD,&d)); 88 PetscCall(VecSetSizes(d,m>n ? n : m,PETSC_DECIDE)); 89 PetscCall(VecSetFromOptions(d)); 90 PetscCall(MatGetDiagonal(C,d)); 91 if (mats_view) { 92 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Diagonal of C:\n")); 93 PetscCall(VecView(d,PETSC_VIEWER_STDOUT_WORLD)); 94 } 95 if (m>n) { 96 PetscCall(MatDiagonalScale(C,NULL,d)); 97 } else { 98 PetscCall(MatDiagonalScale(C,d,NULL)); 99 } 100 if (mats_view) { 101 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Diagonal Scaled C:\n")); 102 PetscCall(MatView(C,PETSC_VIEWER_STDOUT_WORLD)); 103 } 104 105 /* Test MatAXPY(), MatAYPX() and in-place MatConvert() */ 106 PetscCall(MatCreate(PETSC_COMM_WORLD,&B)); 107 PetscCall(MatSetSizes(B,m,n,PETSC_DECIDE,PETSC_DECIDE)); 108 PetscCall(MatSetType(B,MATELEMENTAL)); 109 PetscCall(MatSetFromOptions(B)); 110 PetscCall(MatSetUp(B)); 111 PetscCall(MatGetOwnershipIS(B,&isrows,&iscols)); 112 PetscCall(ISGetLocalSize(isrows,&nrows)); 113 PetscCall(ISGetIndices(isrows,&rows)); 114 PetscCall(ISGetLocalSize(iscols,&ncols)); 115 PetscCall(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 PetscCall(MatSetValues(B,nrows,rows,ncols,cols,v,INSERT_VALUES)); 122 PetscCall(PetscFree(v)); 123 PetscCall(ISRestoreIndices(isrows,&rows)); 124 PetscCall(ISRestoreIndices(iscols,&cols)); 125 PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 126 PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 127 PetscCall(MatAXPY(B,2.5,C,SAME_NONZERO_PATTERN)); 128 PetscCall(MatAYPX(B,3.75,C,SAME_NONZERO_PATTERN)); 129 PetscCall(MatConvert(B,MATDENSE,MAT_INPLACE_MATRIX,&B)); 130 if (mats_view) { 131 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"B after MatAXPY and MatAYPX:\n")); 132 PetscCall(MatView(B,PETSC_VIEWER_STDOUT_WORLD)); 133 } 134 PetscCall(ISDestroy(&isrows)); 135 PetscCall(ISDestroy(&iscols)); 136 PetscCall(MatDestroy(&B)); 137 138 /* Test MatMatTransposeMult(): B = C*C^T */ 139 PetscCall(MatMatTransposeMult(C,C,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&B)); 140 PetscCall(MatScale(C,2.0)); 141 PetscCall(MatMatTransposeMult(C,C,MAT_REUSE_MATRIX,PETSC_DEFAULT,&B)); 142 PetscCall(MatMatTransposeMultEqual(C,C,B,10,&flg)); 143 PetscCheck(flg,PETSC_COMM_WORLD,PETSC_ERR_ARG_INCOMP,"MatMatTransposeMult: B != C*B^T"); 144 145 if (mats_view) { 146 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"C MatMatTransposeMult C:\n")); 147 PetscCall(MatView(B,PETSC_VIEWER_STDOUT_WORLD)); 148 } 149 150 PetscCall(MatDestroy(&Cdense)); 151 PetscCall(PetscFree(v)); 152 PetscCall(MatDestroy(&B)); 153 PetscCall(MatDestroy(&C)); 154 PetscCall(MatDestroy(&Ct)); 155 PetscCall(VecDestroy(&d)); 156 PetscCall(PetscFinalize()); 157 return 0; 158 } 159 160 /*TEST 161 162 test: 163 nsize: 2 164 args: -m 3 -n 2 165 requires: elemental 166 167 test: 168 suffix: 2 169 nsize: 6 170 args: -m 2 -n 3 171 requires: elemental 172 173 test: 174 suffix: 3 175 nsize: 1 176 args: -m 2 -n 3 177 requires: elemental 178 output_file: output/ex39_1.out 179 180 TEST*/ 181