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