1 #include <petscmat.h> 2 3 static char help[] = "Tests MatMatMult with MAT_REUSE_MATRIX and already allocated dense result.\n\n"; 4 5 static PetscErrorCode CheckLocal(Mat A, Mat B, PetscScalar *a, PetscScalar *b) 6 { 7 PetscErrorCode ierr; 8 PetscBool wA = PETSC_FALSE, wB = PETSC_FALSE; 9 10 PetscFunctionBegin; 11 if (a) { 12 const PetscScalar *Aa; 13 ierr = MatDenseGetArrayRead(A,&Aa);CHKERRQ(ierr); 14 wA = (PetscBool)(a != Aa); 15 ierr = MatDenseRestoreArrayRead(A,&Aa);CHKERRQ(ierr); 16 } 17 if (b) { 18 const PetscScalar *Bb; 19 ierr = MatDenseGetArrayRead(B,&Bb);CHKERRQ(ierr); 20 wB = (PetscBool)(b != Bb); 21 ierr = MatDenseRestoreArrayRead(B,&Bb);CHKERRQ(ierr); 22 } 23 if (wA || wB) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Wrong array in A Mat %d, Wrong array in B Mat %d",wA,wB); 24 PetscFunctionReturn(0); 25 } 26 27 int main(int argc,char **args) 28 { 29 Mat X,B,A,T; 30 PetscInt m,n,M = 10,N = 10,K = 5, ldx = 0, ldb = 0; 31 const char *deft = MATAIJ; 32 char mattype[256]; 33 PetscBool flg,symm = PETSC_FALSE,testtt = PETSC_TRUE, testnest = PETSC_TRUE, testtranspose = PETSC_TRUE, testcircular = PETSC_FALSE, local = PETSC_TRUE; 34 PetscScalar *dataX = NULL,*dataB = NULL; 35 PetscScalar *aX,*aB; 36 PetscErrorCode ierr; 37 38 ierr = PetscInitialize(&argc,&args,NULL,help);if (ierr) return ierr; 39 ierr = PetscOptionsGetInt(NULL,NULL,"-N",&N,NULL);CHKERRQ(ierr); 40 ierr = PetscOptionsGetInt(NULL,NULL,"-M",&M,NULL);CHKERRQ(ierr); 41 ierr = PetscOptionsGetInt(NULL,NULL,"-K",&K,NULL);CHKERRQ(ierr); 42 ierr = PetscOptionsGetBool(NULL,NULL,"-symm",&symm,NULL);CHKERRQ(ierr); 43 ierr = PetscOptionsGetBool(NULL,NULL,"-local",&local,NULL);CHKERRQ(ierr); 44 ierr = PetscOptionsGetInt(NULL,NULL,"-ldx",&ldx,NULL);CHKERRQ(ierr); 45 ierr = PetscOptionsGetInt(NULL,NULL,"-ldb",&ldb,NULL);CHKERRQ(ierr); 46 ierr = PetscOptionsGetBool(NULL,NULL,"-testtranspose",&testtranspose,NULL);CHKERRQ(ierr); 47 ierr = PetscOptionsGetBool(NULL,NULL,"-testnest",&testnest,NULL);CHKERRQ(ierr); 48 ierr = PetscOptionsGetBool(NULL,NULL,"-testtt",&testtt,NULL);CHKERRQ(ierr); 49 ierr = PetscOptionsGetBool(NULL,NULL,"-testcircular",&testcircular,NULL);CHKERRQ(ierr); 50 #if defined(PETSC_USE_COMPLEX) 51 testtranspose = PETSC_FALSE; 52 testtt = PETSC_FALSE; 53 #endif 54 ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr); 55 ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,M,N);CHKERRQ(ierr); 56 ierr = MatSetType(A,MATAIJ);CHKERRQ(ierr); 57 ierr = MatSetUp(A);CHKERRQ(ierr); 58 ierr = MatSetRandom(A,NULL);CHKERRQ(ierr); 59 if (M==N && symm) { 60 Mat AT; 61 62 ierr = MatHermitianTranspose(A,MAT_INITIAL_MATRIX,&AT);CHKERRQ(ierr); 63 ierr = MatAXPY(A,1.0,AT,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr); 64 ierr = MatDestroy(&AT);CHKERRQ(ierr); 65 ierr = MatSetOption(A,MAT_HERMITIAN,PETSC_TRUE);CHKERRQ(ierr); 66 } 67 ierr = MatViewFromOptions(A,NULL,"-A_init_view");CHKERRQ(ierr); 68 ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"","","");CHKERRQ(ierr); 69 ierr = PetscOptionsFList("-mat_type","Matrix type","MatSetType",MatList,deft,mattype,256,&flg);CHKERRQ(ierr); 70 ierr = PetscOptionsEnd();CHKERRQ(ierr); 71 if (flg) { 72 Mat A2; 73 74 ierr = MatConvert(A,mattype,MAT_INITIAL_MATRIX,&A2);CHKERRQ(ierr); 75 ierr = MatMultEqual(A,A2,10,&flg);CHKERRQ(ierr); 76 if (!flg) { 77 Mat AE,A2E; 78 79 ierr = PetscPrintf(PETSC_COMM_WORLD,"Error with convert\n");CHKERRQ(ierr); 80 ierr = MatComputeOperator(A,MATDENSE,&AE);CHKERRQ(ierr); 81 ierr = MatComputeOperator(A2,MATDENSE,&A2E);CHKERRQ(ierr); 82 ierr = MatView(AE,NULL);CHKERRQ(ierr); 83 ierr = MatView(A2E,NULL);CHKERRQ(ierr); 84 ierr = MatAXPY(A2E,-1.0,A,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 85 ierr = MatView(A2E,NULL);CHKERRQ(ierr); 86 ierr = MatDestroy(&A2E);CHKERRQ(ierr); 87 ierr = MatDestroy(&AE);CHKERRQ(ierr); 88 } 89 ierr = MatDestroy(&A);CHKERRQ(ierr); 90 A = A2; 91 } 92 ierr = MatViewFromOptions(A,NULL,"-A_view");CHKERRQ(ierr); 93 94 ierr = MatGetLocalSize(A,&m,&n);CHKERRQ(ierr); 95 if (local) { 96 ierr = PetscMalloc1(PetscMax((m+ldx)*K,1),&dataX);CHKERRQ(ierr); 97 ierr = PetscMalloc1(PetscMax((n+ldb)*K,1),&dataB);CHKERRQ(ierr); 98 } 99 ierr = MatCreateDense(PETSC_COMM_WORLD,n,PETSC_DECIDE,N,K,dataB,&B);CHKERRQ(ierr); 100 ierr = MatCreateDense(PETSC_COMM_WORLD,m,PETSC_DECIDE,M,K,dataX,&X);CHKERRQ(ierr); 101 102 /* store pointer to dense data for testing */ 103 ierr = MatDenseGetArrayRead(B,(const PetscScalar**)&dataB);CHKERRQ(ierr); 104 ierr = MatDenseGetArrayRead(X,(const PetscScalar**)&dataX);CHKERRQ(ierr); 105 aX = dataX; 106 aB = dataB; 107 ierr = MatDenseRestoreArrayRead(B,(const PetscScalar**)&dataB);CHKERRQ(ierr); 108 ierr = MatDenseRestoreArrayRead(X,(const PetscScalar**)&dataX);CHKERRQ(ierr); 109 if (local) { 110 dataX = aX; 111 dataB = aB; 112 } 113 ierr = MatSetRandom(B,NULL);CHKERRQ(ierr); 114 115 /* Test reusing a previously allocated dense buffer */ 116 ierr = MatMatMult(A,B,MAT_REUSE_MATRIX,PETSC_DEFAULT,&X);CHKERRQ(ierr); 117 ierr = CheckLocal(B,X,aB,aX);CHKERRQ(ierr); 118 ierr = MatMatMultEqual(A,B,X,10,&flg);CHKERRQ(ierr); 119 if (!flg) { 120 ierr = PetscPrintf(PETSC_COMM_WORLD,"Error with reusage\n");CHKERRQ(ierr); 121 ierr = MatMatMult(A,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&T);CHKERRQ(ierr); 122 ierr = MatAXPY(T,-1.0,X,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 123 ierr = MatView(T,NULL);CHKERRQ(ierr); 124 ierr = MatDestroy(&T);CHKERRQ(ierr); 125 } 126 127 if (testnest) { /* test with MatNest */ 128 Mat NA; 129 130 ierr = MatCreateNest(PETSC_COMM_WORLD,1,NULL,1,NULL,&A,&NA);CHKERRQ(ierr); 131 ierr = MatViewFromOptions(NA,NULL,"-NA_view");CHKERRQ(ierr); 132 ierr = MatMatMult(NA,B,MAT_REUSE_MATRIX,PETSC_DEFAULT,&X);CHKERRQ(ierr); 133 ierr = CheckLocal(B,X,aB,aX);CHKERRQ(ierr); 134 ierr = MatMatMultEqual(NA,B,X,10,&flg);CHKERRQ(ierr); 135 if (!flg) { 136 ierr = PetscPrintf(PETSC_COMM_WORLD,"Error with Nest\n");CHKERRQ(ierr); 137 ierr = MatMatMult(NA,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&T);CHKERRQ(ierr); 138 ierr = MatAXPY(T,-1.0,X,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 139 ierr = MatView(T,NULL);CHKERRQ(ierr); 140 ierr = MatDestroy(&T);CHKERRQ(ierr); 141 } 142 ierr = MatDestroy(&NA);CHKERRQ(ierr); 143 } 144 145 if (testtranspose) { /* test with Transpose */ 146 Mat TA; 147 148 ierr = MatCreateHermitianTranspose(A,&TA);CHKERRQ(ierr); 149 ierr = MatMatMult(TA,X,MAT_REUSE_MATRIX,PETSC_DEFAULT,&B);CHKERRQ(ierr); 150 ierr = CheckLocal(B,X,aB,aX);CHKERRQ(ierr); 151 ierr = MatMatMultEqual(TA,X,B,10,&flg);CHKERRQ(ierr); 152 if (!flg) { 153 ierr = PetscPrintf(PETSC_COMM_WORLD,"Error with Transpose\n");CHKERRQ(ierr); 154 ierr = MatMatMult(TA,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&T);CHKERRQ(ierr); 155 ierr = MatAXPY(T,-1.0,X,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 156 ierr = MatView(T,NULL);CHKERRQ(ierr); 157 ierr = MatDestroy(&T);CHKERRQ(ierr); 158 } 159 ierr = MatDestroy(&TA);CHKERRQ(ierr); 160 } 161 162 if (testtt) { /* test with Transpose(Transpose) */ 163 Mat TA, TTA; 164 165 ierr = MatCreateHermitianTranspose(A,&TA);CHKERRQ(ierr); 166 ierr = MatCreateHermitianTranspose(TA,&TTA);CHKERRQ(ierr); 167 ierr = MatMatMult(TTA,B,MAT_REUSE_MATRIX,PETSC_DEFAULT,&X);CHKERRQ(ierr); 168 ierr = CheckLocal(B,X,aB,aX);CHKERRQ(ierr); 169 ierr = MatMatMultEqual(TTA,B,X,10,&flg);CHKERRQ(ierr); 170 if (!flg) { 171 ierr = PetscPrintf(PETSC_COMM_WORLD,"Error with Transpose(Transpose)\n");CHKERRQ(ierr); 172 ierr = MatMatMult(TTA,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&T);CHKERRQ(ierr); 173 ierr = MatAXPY(T,-1.0,X,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 174 ierr = MatView(T,NULL);CHKERRQ(ierr); 175 ierr = MatDestroy(&T);CHKERRQ(ierr); 176 } 177 ierr = MatDestroy(&TA);CHKERRQ(ierr); 178 ierr = MatDestroy(&TTA);CHKERRQ(ierr); 179 } 180 181 if (testcircular) { /* test circular */ 182 Mat AB; 183 184 ierr = MatMatMult(A,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&AB);CHKERRQ(ierr); 185 ierr = MatMatMult(A,B,MAT_REUSE_MATRIX,PETSC_DEFAULT,&X);CHKERRQ(ierr); 186 ierr = CheckLocal(B,X,aB,aX);CHKERRQ(ierr); 187 if (M == N && N == K) { 188 ierr = MatMatMult(A,X,MAT_REUSE_MATRIX,PETSC_DEFAULT,&B);CHKERRQ(ierr); 189 } else { 190 ierr = MatTransposeMatMult(A,X,MAT_REUSE_MATRIX,PETSC_DEFAULT,&B);CHKERRQ(ierr); 191 } 192 ierr = CheckLocal(B,X,aB,aX);CHKERRQ(ierr); 193 ierr = MatDestroy(&AB);CHKERRQ(ierr); 194 } 195 ierr = MatDestroy(&X);CHKERRQ(ierr); 196 ierr = MatDestroy(&B);CHKERRQ(ierr); 197 ierr = MatDestroy(&A);CHKERRQ(ierr); 198 ierr = PetscFree(dataX);CHKERRQ(ierr); 199 ierr = PetscFree(dataB);CHKERRQ(ierr); 200 ierr = PetscFinalize(); 201 return ierr; 202 } 203 204 /*TEST 205 206 test: 207 suffix: 1 208 args: -local {{0 1}} 209 210 test: 211 output_file: output/ex70_1.out 212 nsize: 2 213 suffix: 1_par 214 args: -testtranspose 0 -local {{0 1}} 215 216 test: 217 TODO: MPIAIJ x MPIDENSE broken for MatTransposeMatMult 218 output_file: output/ex70_1.out 219 nsize: 2 220 suffix: 1_par_broken 221 args: -testtranspose -local {{0 1}} 222 223 test: 224 output_file: output/ex70_1.out 225 suffix: 2 226 nsize: 1 227 args: -M {{7 11}} -N {{12 9}} -K {{1 3}} -local {{0 1}} 228 229 test: 230 output_file: output/ex70_1.out 231 suffix: 2_par 232 nsize: 2 233 args: -M {{7 11}} -N {{12 9}} -K {{1 3}} -local {{0 1}} -testcircular 0 234 235 test: 236 TODO: MatTransposeMatMultSymbolic_SeqAIJ_SeqDense plays with the destroy routine 237 output_file: output/ex70_1.out 238 suffix: 2_broken 239 nsize: 2 240 args: -M {{7 11}} -N {{12 9}} -K {{1 3}} -local {{0 1}} -testcircular 1 241 242 test: 243 output_file: output/ex70_1.out 244 suffix: 3 245 nsize: 1 246 args: -M 13 -N 13 -K {{1 3}} -local {{0 1}} -mat_type sbaij -symm -testtranspose 0 247 248 test: 249 output_file: output/ex70_1.out 250 suffix: 4 251 nsize: 1 252 args: -M 3 -N 3 -K 3 -local {{0 1}} -testcircular 253 254 test: 255 output_file: output/ex70_1.out 256 suffix: 5 257 nsize: {{2 4}} 258 args: -M 3 -N 3 -K 3 -local 1 -testcircular -testtranspose 0 259 260 test: 261 TODO: MatCreateDense broken with some processors not having local rows 262 output_file: output/ex70_1.out 263 suffix: 5_broken 264 nsize: {{2 4}} 265 args: -M 3 -N 3 -K 3 -local 0 -testcircular -testtranspose 0 266 267 test: 268 output_file: output/ex70_1.out 269 suffix: 6 270 nsize: 1 271 args: -M {{1 3}} -N {{2 5}} -K {{1 2}} -local {{0 1}} -testcircular -testtranspose 0 272 273 test: 274 TODO: MatTransposeMatMultSymbolic_SeqAIJ_SeqDense plays with the destroy routine 275 output_file: output/ex70_1.out 276 suffix: 6_broken 277 nsize: 1 278 args: -M {{1 3}} -N {{2 5}} -K {{1 2}} -local {{0 1}} -testcircular -testtranspose 279 280 test: 281 output_file: output/ex70_1.out 282 suffix: 7 283 nsize: 1 284 args: -M 13 -N 13 -K {{1 3}} -local {{0 1}} -mat_type dense -testnest 0 -testcircular 285 286 test: 287 TODO: NEST x DENSE with dense nested matrices seems broken in this case 288 output_file: output/ex70_1.out 289 suffix: 7_broken 290 nsize: 1 291 args: -M 13 -N 13 -K {{1 3}} -local {{0 1}} -mat_type dense -testnest -testcircular 292 TEST*/ 293