1 #include <petscmat.h> 2 3 static char help[] = "Tests MatMat operations with MAT_REUSE_MATRIX and already allocated dense result.\n\n"; 4 5 static PetscScalar MAGIC_NUMBER = 12345; 6 7 static PetscErrorCode CheckLocal(Mat A, Mat B, PetscScalar *a, PetscScalar *b) 8 { 9 PetscErrorCode ierr; 10 PetscBool wA = PETSC_FALSE, wB = PETSC_FALSE; 11 PetscBool wAv = PETSC_FALSE, wBv = PETSC_FALSE; 12 PetscInt lda,i,j,m,n; 13 14 PetscFunctionBegin; 15 if (a) { 16 const PetscScalar *Aa; 17 ierr = MatDenseGetArrayRead(A,&Aa);CHKERRQ(ierr); 18 wA = (PetscBool)(a != Aa); 19 ierr = MatDenseGetLDA(A,&lda);CHKERRQ(ierr); 20 ierr = MatGetLocalSize(A,&m,&n);CHKERRQ(ierr); 21 for (j=0;j<n;j++) { 22 for (i=m;i<lda;i++) { 23 if (Aa[j*lda +i] != MAGIC_NUMBER) wAv = PETSC_TRUE; 24 } 25 } 26 ierr = MatDenseRestoreArrayRead(A,&Aa);CHKERRQ(ierr); 27 } 28 if (b) { 29 const PetscScalar *Bb; 30 ierr = MatDenseGetArrayRead(B,&Bb);CHKERRQ(ierr); 31 wB = (PetscBool)(b != Bb); 32 ierr = MatDenseGetLDA(B,&lda);CHKERRQ(ierr); 33 ierr = MatGetLocalSize(B,&m,&n);CHKERRQ(ierr); 34 for (j=0;j<n;j++) { 35 for (i=m;i<lda;i++) { 36 if (Bb[j*lda +i] != MAGIC_NUMBER) wBv = PETSC_TRUE; 37 } 38 } 39 ierr = MatDenseRestoreArrayRead(B,&Bb);CHKERRQ(ierr); 40 } 41 if (wA || wB) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Wrong array in first Mat? %d, Wrong array in second Mat? %d",wA,wB); 42 if (wAv || wBv) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Wrong data in first Mat? %d, Wrong data in second Mat? %d",wAv,wBv); 43 PetscFunctionReturn(0); 44 } 45 46 typedef struct { 47 Mat A; 48 Mat P; 49 Mat R; 50 } proj_data; 51 52 PetscErrorCode proj_destroy(void *ctx) 53 { 54 proj_data *userdata = (proj_data*)ctx; 55 PetscErrorCode ierr; 56 57 PetscFunctionBegin; 58 if (!userdata) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Missing userdata"); 59 ierr = MatDestroy(&userdata->A);CHKERRQ(ierr); 60 ierr = MatDestroy(&userdata->P);CHKERRQ(ierr); 61 ierr = MatDestroy(&userdata->R);CHKERRQ(ierr); 62 ierr = PetscFree(userdata);CHKERRQ(ierr); 63 PetscFunctionReturn(0); 64 } 65 66 PetscErrorCode proj_mult(Mat S, Vec X, Vec Y) 67 { 68 Mat A,R,P; 69 Vec Ax,Ay; 70 Vec Px,Py; 71 proj_data *userdata; 72 PetscErrorCode ierr; 73 74 PetscFunctionBegin; 75 ierr = MatShellGetContext(S,(void**)&userdata);CHKERRQ(ierr); 76 if (!userdata) SETERRQ(PetscObjectComm((PetscObject)S),PETSC_ERR_PLIB,"Missing userdata"); 77 A = userdata->A; 78 R = userdata->R; 79 P = userdata->P; 80 if (!A) SETERRQ(PetscObjectComm((PetscObject)S),PETSC_ERR_PLIB,"Missing matrix"); 81 if (!R && !P) SETERRQ(PetscObjectComm((PetscObject)S),PETSC_ERR_PLIB,"Missing projectors"); 82 if (R && P) SETERRQ(PetscObjectComm((PetscObject)S),PETSC_ERR_PLIB,"Both projectors"); 83 ierr = MatCreateVecs(A,&Ax,&Ay);CHKERRQ(ierr); 84 if (R) { 85 ierr = MatCreateVecs(R,&Py,&Px);CHKERRQ(ierr); 86 } else { 87 ierr = MatCreateVecs(P,&Px,&Py);CHKERRQ(ierr); 88 } 89 ierr = VecCopy(X,Px);CHKERRQ(ierr); 90 if (P) { 91 ierr = MatMult(P,Px,Py);CHKERRQ(ierr); 92 } else { 93 ierr = MatMultTranspose(R,Px,Py);CHKERRQ(ierr); 94 } 95 ierr = VecCopy(Py,Ax);CHKERRQ(ierr); 96 ierr = MatMult(A,Ax,Ay);CHKERRQ(ierr); 97 ierr = VecCopy(Ay,Py);CHKERRQ(ierr); 98 if (P) { 99 ierr = MatMultTranspose(P,Py,Px);CHKERRQ(ierr); 100 } else { 101 ierr = MatMult(R,Py,Px);CHKERRQ(ierr); 102 } 103 ierr = VecCopy(Px,Y);CHKERRQ(ierr); 104 ierr = VecDestroy(&Px);CHKERRQ(ierr); 105 ierr = VecDestroy(&Py);CHKERRQ(ierr); 106 ierr = VecDestroy(&Ax);CHKERRQ(ierr); 107 ierr = VecDestroy(&Ay);CHKERRQ(ierr); 108 PetscFunctionReturn(0); 109 } 110 111 PetscErrorCode MyPtShellPMultSymbolic(Mat S, Mat P, Mat PtAP, void** ctx) 112 { 113 PetscErrorCode ierr; 114 proj_data *userdata; 115 116 PetscFunctionBegin; 117 ierr = PetscNew(&userdata);CHKERRQ(ierr); 118 ierr = MatShellSetContext(PtAP,(void*)userdata);CHKERRQ(ierr); 119 *ctx = (void *)userdata; 120 PetscFunctionReturn(0); 121 } 122 123 PetscErrorCode MyPtShellPMultNumeric(Mat S, Mat P, Mat PtAP, void *ctx) 124 { 125 Mat A; 126 PetscErrorCode ierr; 127 proj_data *userdata = (proj_data*)ctx; 128 129 PetscFunctionBegin; 130 ierr = MatShellGetContext(S,(void**)&A);CHKERRQ(ierr); 131 ierr = PetscObjectReference((PetscObject)A);CHKERRQ(ierr); 132 ierr = PetscObjectReference((PetscObject)P);CHKERRQ(ierr); 133 ierr = MatDestroy(&userdata->A);CHKERRQ(ierr); 134 ierr = MatDestroy(&userdata->P);CHKERRQ(ierr); 135 ierr = MatDestroy(&userdata->R);CHKERRQ(ierr); 136 userdata->A = A; 137 userdata->P = P; 138 ierr = MatShellSetOperation(PtAP,MATOP_MULT,(void (*)(void))proj_mult);CHKERRQ(ierr); 139 ierr = MatSetUp(PtAP);CHKERRQ(ierr); 140 ierr = MatAssemblyBegin(PtAP,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 141 ierr = MatAssemblyEnd(PtAP,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 142 PetscFunctionReturn(0); 143 } 144 145 PetscErrorCode MyRShellRtMultSymbolic(Mat S, Mat R, Mat RARt, void **ctx) 146 { 147 PetscErrorCode ierr; 148 proj_data *userdata; 149 150 PetscFunctionBegin; 151 ierr = PetscNew(&userdata);CHKERRQ(ierr); 152 ierr = MatShellSetContext(RARt,(void*)userdata);CHKERRQ(ierr); 153 *ctx = (void *)userdata; 154 PetscFunctionReturn(0); 155 } 156 157 PetscErrorCode MyRShellRtMultNumeric(Mat S, Mat R, Mat RARt, void *ctx) 158 { 159 Mat A; 160 PetscErrorCode ierr; 161 proj_data *userdata = (proj_data*)ctx; 162 163 PetscFunctionBegin; 164 ierr = MatShellGetContext(S,(void**)&A);CHKERRQ(ierr); 165 ierr = PetscObjectReference((PetscObject)A);CHKERRQ(ierr); 166 ierr = PetscObjectReference((PetscObject)R);CHKERRQ(ierr); 167 ierr = MatDestroy(&userdata->A);CHKERRQ(ierr); 168 ierr = MatDestroy(&userdata->P);CHKERRQ(ierr); 169 ierr = MatDestroy(&userdata->R);CHKERRQ(ierr); 170 userdata->A = A; 171 userdata->R = R; 172 ierr = MatShellSetOperation(RARt,MATOP_MULT,(void (*)(void))proj_mult);CHKERRQ(ierr); 173 ierr = MatSetUp(RARt);CHKERRQ(ierr); 174 ierr = MatAssemblyBegin(RARt,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 175 ierr = MatAssemblyEnd(RARt,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 176 PetscFunctionReturn(0); 177 } 178 179 PetscErrorCode MyMatShellMatMultNumeric(Mat S, Mat B, Mat C, void *ctx) 180 { 181 PetscErrorCode ierr; 182 Mat A; 183 184 PetscFunctionBegin; 185 ierr = MatShellGetContext(S,(void**)&A);CHKERRQ(ierr); 186 ierr = MatMatMult(A,B,MAT_REUSE_MATRIX,PETSC_DEFAULT,&C);CHKERRQ(ierr); 187 PetscFunctionReturn(0); 188 } 189 190 PetscErrorCode MyMatTransposeShellMatMultNumeric(Mat S, Mat B, Mat C, void *ctx) 191 { 192 PetscErrorCode ierr; 193 Mat A; 194 195 PetscFunctionBegin; 196 ierr = MatShellGetContext(S,(void**)&A);CHKERRQ(ierr); 197 ierr = MatTransposeMatMult(A,B,MAT_REUSE_MATRIX,PETSC_DEFAULT,&C);CHKERRQ(ierr); 198 PetscFunctionReturn(0); 199 } 200 201 PetscErrorCode MyMatShellMatTransposeMultNumeric(Mat S, Mat B, Mat C, void *ctx) 202 { 203 PetscErrorCode ierr; 204 Mat A; 205 206 PetscFunctionBegin; 207 ierr = MatShellGetContext(S,(void**)&A);CHKERRQ(ierr); 208 ierr = MatMatTransposeMult(A,B,MAT_REUSE_MATRIX,PETSC_DEFAULT,&C);CHKERRQ(ierr); 209 PetscFunctionReturn(0); 210 } 211 212 int main(int argc,char **args) 213 { 214 Mat X,B,A,Bt,T,T2,PtAP = NULL,RARt = NULL, R = NULL; 215 Vec r,l,rs,ls; 216 PetscInt m,n,k,M = 10,N = 10,K = 5, ldx = 3, ldb = 5, ldr = 4; 217 const char *deft = MATAIJ; 218 char mattype[256]; 219 PetscBool flg,symm = PETSC_FALSE,testtt = PETSC_TRUE, testnest = PETSC_TRUE, testtranspose = PETSC_TRUE, testcircular = PETSC_FALSE, local = PETSC_TRUE; 220 PetscBool testhtranspose = PETSC_TRUE; 221 PetscBool xgpu = PETSC_FALSE, bgpu = PETSC_FALSE, testshellops = PETSC_FALSE, testproj = PETSC_TRUE, testrart = PETSC_TRUE, testmatmatt = PETSC_TRUE, testmattmat = PETSC_TRUE; 222 PetscScalar *dataX = NULL,*dataB = NULL, *dataR = NULL, *dataBt = NULL; 223 PetscScalar *aX,*aB,*aBt; 224 PetscReal err; 225 PetscErrorCode ierr; 226 227 ierr = PetscInitialize(&argc,&args,NULL,help);if (ierr) return ierr; 228 ierr = PetscOptionsGetInt(NULL,NULL,"-N",&N,NULL);CHKERRQ(ierr); 229 ierr = PetscOptionsGetInt(NULL,NULL,"-M",&M,NULL);CHKERRQ(ierr); 230 ierr = PetscOptionsGetInt(NULL,NULL,"-K",&K,NULL);CHKERRQ(ierr); 231 ierr = PetscOptionsGetBool(NULL,NULL,"-symm",&symm,NULL);CHKERRQ(ierr); 232 ierr = PetscOptionsGetBool(NULL,NULL,"-local",&local,NULL);CHKERRQ(ierr); 233 ierr = PetscOptionsGetInt(NULL,NULL,"-ldx",&ldx,NULL);CHKERRQ(ierr); 234 ierr = PetscOptionsGetInt(NULL,NULL,"-ldb",&ldb,NULL);CHKERRQ(ierr); 235 ierr = PetscOptionsGetInt(NULL,NULL,"-ldr",&ldr,NULL);CHKERRQ(ierr); 236 ierr = PetscOptionsGetBool(NULL,NULL,"-testtranspose",&testtranspose,NULL);CHKERRQ(ierr); 237 ierr = PetscOptionsGetBool(NULL,NULL,"-testnest",&testnest,NULL);CHKERRQ(ierr); 238 ierr = PetscOptionsGetBool(NULL,NULL,"-testtt",&testtt,NULL);CHKERRQ(ierr); 239 ierr = PetscOptionsGetBool(NULL,NULL,"-testcircular",&testcircular,NULL);CHKERRQ(ierr); 240 ierr = PetscOptionsGetBool(NULL,NULL,"-testshellops",&testshellops,NULL);CHKERRQ(ierr); 241 ierr = PetscOptionsGetBool(NULL,NULL,"-testproj",&testproj,NULL);CHKERRQ(ierr); 242 ierr = PetscOptionsGetBool(NULL,NULL,"-testrart",&testrart,NULL);CHKERRQ(ierr); 243 ierr = PetscOptionsGetBool(NULL,NULL,"-testmatmatt",&testmatmatt,NULL);CHKERRQ(ierr); 244 ierr = PetscOptionsGetBool(NULL,NULL,"-testmattmat",&testmattmat,NULL);CHKERRQ(ierr); 245 ierr = PetscOptionsGetBool(NULL,NULL,"-xgpu",&xgpu,NULL);CHKERRQ(ierr); 246 ierr = PetscOptionsGetBool(NULL,NULL,"-bgpu",&bgpu,NULL);CHKERRQ(ierr); 247 ierr = PetscOptionsGetScalar(NULL,NULL,"-magic_number",&MAGIC_NUMBER,NULL);CHKERRQ(ierr); 248 if (M != N) testproj = PETSC_FALSE; 249 250 ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr); 251 ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,M,N);CHKERRQ(ierr); 252 ierr = MatSetType(A,MATAIJ);CHKERRQ(ierr); 253 ierr = MatSetUp(A);CHKERRQ(ierr); 254 ierr = MatSetRandom(A,NULL);CHKERRQ(ierr); 255 if (M==N && symm) { 256 Mat AT; 257 258 ierr = MatTranspose(A,MAT_INITIAL_MATRIX,&AT);CHKERRQ(ierr); 259 ierr = MatAXPY(A,1.0,AT,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr); 260 ierr = MatDestroy(&AT);CHKERRQ(ierr); 261 ierr = MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);CHKERRQ(ierr); 262 } 263 ierr = MatViewFromOptions(A,NULL,"-A_init_view");CHKERRQ(ierr); 264 ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"","","");CHKERRQ(ierr); 265 ierr = PetscOptionsFList("-A_mat_type","Matrix type","MatSetType",MatList,deft,mattype,256,&flg);CHKERRQ(ierr); 266 ierr = PetscOptionsEnd();CHKERRQ(ierr); 267 if (flg) { 268 Mat A2; 269 270 ierr = MatDuplicate(A,MAT_COPY_VALUES,&A2);CHKERRQ(ierr); 271 ierr = MatConvert(A,mattype,MAT_INPLACE_MATRIX,&A);CHKERRQ(ierr); 272 ierr = MatMultEqual(A,A2,10,&flg);CHKERRQ(ierr); 273 if (!flg) { 274 Mat AE,A2E; 275 276 ierr = PetscPrintf(PETSC_COMM_WORLD,"Error with convert\n");CHKERRQ(ierr); 277 ierr = MatComputeOperator(A,MATDENSE,&AE);CHKERRQ(ierr); 278 ierr = MatComputeOperator(A2,MATDENSE,&A2E);CHKERRQ(ierr); 279 ierr = MatView(AE,NULL);CHKERRQ(ierr); 280 ierr = MatView(A2E,NULL);CHKERRQ(ierr); 281 ierr = MatAXPY(A2E,-1.0,A,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 282 ierr = MatView(A2E,NULL);CHKERRQ(ierr); 283 ierr = MatDestroy(&A2E);CHKERRQ(ierr); 284 ierr = MatDestroy(&AE);CHKERRQ(ierr); 285 } 286 ierr = MatDestroy(&A2);CHKERRQ(ierr); 287 } 288 ierr = MatViewFromOptions(A,NULL,"-A_view");CHKERRQ(ierr); 289 290 ierr = MatGetLocalSize(A,&m,&n);CHKERRQ(ierr); 291 if (local) { 292 PetscInt i; 293 294 ierr = PetscMalloc1((m+ldx)*K,&dataX);CHKERRQ(ierr); 295 ierr = PetscMalloc1((n+ldb)*K,&dataB);CHKERRQ(ierr); 296 for (i=0;i<(m+ldx)*K;i++) dataX[i] = MAGIC_NUMBER; 297 for (i=0;i<(n+ldb)*K;i++) dataB[i] = MAGIC_NUMBER; 298 } 299 ierr = MatCreateDense(PETSC_COMM_WORLD,n,PETSC_DECIDE,N,K,dataB,&B);CHKERRQ(ierr); 300 ierr = MatCreateDense(PETSC_COMM_WORLD,m,PETSC_DECIDE,M,K,dataX,&X);CHKERRQ(ierr); 301 if (local) { 302 ierr = MatDenseSetLDA(X,m+ldx);CHKERRQ(ierr); 303 ierr = MatDenseSetLDA(B,n+ldb);CHKERRQ(ierr); 304 } 305 ierr = MatGetLocalSize(B,NULL,&k);CHKERRQ(ierr); 306 if (local) { 307 PetscInt i; 308 309 ierr = PetscMalloc1((k+ldr)*N,&dataBt);CHKERRQ(ierr); 310 for (i=0;i<(k+ldr)*N;i++) dataBt[i] = MAGIC_NUMBER; 311 } 312 ierr = MatCreateDense(PETSC_COMM_WORLD,k,n,K,N,dataBt,&Bt);CHKERRQ(ierr); 313 if (local) { 314 ierr = MatDenseSetLDA(Bt,k+ldr);CHKERRQ(ierr); 315 } 316 317 /* store pointer to dense data for testing */ 318 ierr = MatDenseGetArrayRead(B,(const PetscScalar**)&dataB);CHKERRQ(ierr); 319 ierr = MatDenseGetArrayRead(X,(const PetscScalar**)&dataX);CHKERRQ(ierr); 320 ierr = MatDenseGetArrayRead(Bt,(const PetscScalar**)&dataBt);CHKERRQ(ierr); 321 aX = dataX; 322 aB = dataB; 323 aBt = dataBt; 324 ierr = MatDenseRestoreArrayRead(Bt,(const PetscScalar**)&dataBt);CHKERRQ(ierr); 325 ierr = MatDenseRestoreArrayRead(B,(const PetscScalar**)&dataB);CHKERRQ(ierr); 326 ierr = MatDenseRestoreArrayRead(X,(const PetscScalar**)&dataX);CHKERRQ(ierr); 327 if (local) { 328 dataX = aX; 329 dataB = aB; 330 dataBt = aBt; 331 } 332 333 ierr = MatSetRandom(X,NULL);CHKERRQ(ierr); 334 ierr = MatSetRandom(B,NULL);CHKERRQ(ierr); 335 ierr = MatSetRandom(Bt,NULL);CHKERRQ(ierr); 336 ierr = CheckLocal(X,NULL,aX,NULL);CHKERRQ(ierr); 337 ierr = CheckLocal(Bt,B,aBt,aB);CHKERRQ(ierr); 338 339 /* convert to CUDA if needed */ 340 if (bgpu) { 341 ierr = MatConvert(B,MATDENSECUDA,MAT_INPLACE_MATRIX,&B);CHKERRQ(ierr); 342 ierr = MatConvert(Bt,MATDENSECUDA,MAT_INPLACE_MATRIX,&Bt);CHKERRQ(ierr); 343 } 344 if (xgpu) { 345 ierr = MatConvert(X,MATDENSECUDA,MAT_INPLACE_MATRIX,&X);CHKERRQ(ierr); 346 } 347 ierr = CheckLocal(B,X,aB,aX);CHKERRQ(ierr); 348 349 /* Test MatDenseGetSubMatrix */ 350 { 351 Mat B2,T3,T4; 352 353 ierr = MatDuplicate(B,MAT_COPY_VALUES,&B2);CHKERRQ(ierr); 354 ierr = MatDuplicate(B,MAT_DO_NOT_COPY_VALUES,&T4);CHKERRQ(ierr); 355 ierr = MatSetRandom(T4,NULL);CHKERRQ(ierr); 356 ierr = MatAXPY(B2,1.0,T4,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 357 ierr = MatDenseGetSubMatrix(B,PetscMin(1,K),PetscMin(2,K),&T);CHKERRQ(ierr); 358 ierr = MatDenseGetSubMatrix(T4,PetscMin(1,K),PetscMin(2,K),&T2);CHKERRQ(ierr); 359 ierr = MatDenseGetSubMatrix(B2,PetscMin(1,K),PetscMin(2,K),&T3);CHKERRQ(ierr); 360 ierr = MatAXPY(T,1.0,T2,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 361 ierr = MatAXPY(T3,-1.0,T,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 362 ierr = MatNorm(T3,NORM_FROBENIUS,&err);CHKERRQ(ierr); 363 if (err > PETSC_SMALL) { 364 ierr = PetscPrintf(PETSC_COMM_WORLD,"Error with MatDenseGetSubMatrix\n");CHKERRQ(ierr); 365 ierr = MatView(T3,NULL);CHKERRQ(ierr); 366 } 367 ierr = MatDenseRestoreSubMatrix(B,&T);CHKERRQ(ierr); 368 ierr = MatDenseRestoreSubMatrix(T4,&T2);CHKERRQ(ierr); 369 ierr = MatDenseRestoreSubMatrix(B2,&T3);CHKERRQ(ierr); 370 ierr = CheckLocal(B,NULL,aB,NULL);CHKERRQ(ierr); 371 ierr = MatDestroy(&B2);CHKERRQ(ierr); 372 ierr = MatDestroy(&T4);CHKERRQ(ierr); 373 } 374 375 /* Test reusing a previously allocated dense buffer */ 376 ierr = MatMatMult(A,B,MAT_REUSE_MATRIX,PETSC_DEFAULT,&X);CHKERRQ(ierr); 377 ierr = CheckLocal(B,X,aB,aX);CHKERRQ(ierr); 378 ierr = MatMatMultEqual(A,B,X,10,&flg);CHKERRQ(ierr); 379 if (!flg) { 380 ierr = PetscPrintf(PETSC_COMM_WORLD,"Error with reusage\n");CHKERRQ(ierr); 381 ierr = MatMatMult(A,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&T);CHKERRQ(ierr); 382 ierr = MatAXPY(T,-1.0,X,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 383 ierr = MatView(T,NULL);CHKERRQ(ierr); 384 ierr = MatDestroy(&T);CHKERRQ(ierr); 385 } 386 387 /* Test MatTransposeMat and MatMatTranspose */ 388 if (testmattmat) { 389 ierr = MatTransposeMatMult(A,X,MAT_REUSE_MATRIX,PETSC_DEFAULT,&B);CHKERRQ(ierr); 390 ierr = CheckLocal(B,X,aB,aX);CHKERRQ(ierr); 391 ierr = MatTransposeMatMultEqual(A,X,B,10,&flg);CHKERRQ(ierr); 392 if (!flg) { 393 ierr = PetscPrintf(PETSC_COMM_WORLD,"Error with reusage (MatTransposeMat)\n");CHKERRQ(ierr); 394 ierr = MatTransposeMatMult(A,X,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&B);CHKERRQ(ierr); 395 ierr = MatAXPY(T,-1.0,B,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 396 ierr = MatView(T,NULL);CHKERRQ(ierr); 397 ierr = MatDestroy(&T);CHKERRQ(ierr); 398 } 399 } 400 if (testmatmatt) { 401 ierr = MatMatTransposeMult(A,Bt,MAT_REUSE_MATRIX,PETSC_DEFAULT,&X);CHKERRQ(ierr); 402 ierr = CheckLocal(Bt,X,aBt,aX);CHKERRQ(ierr); 403 ierr = MatMatTransposeMultEqual(A,Bt,X,10,&flg);CHKERRQ(ierr); 404 if (!flg) { 405 ierr = PetscPrintf(PETSC_COMM_WORLD,"Error with reusage (MatMatTranspose)\n");CHKERRQ(ierr); 406 ierr = MatMatTransposeMult(A,Bt,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&T);CHKERRQ(ierr); 407 ierr = MatAXPY(T,-1.0,X,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 408 ierr = MatView(T,NULL);CHKERRQ(ierr); 409 ierr = MatDestroy(&T);CHKERRQ(ierr); 410 } 411 } 412 413 /* Test projection operations (PtAP and RARt) */ 414 if (testproj) { 415 ierr = MatPtAP(A,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&PtAP);CHKERRQ(ierr); 416 ierr = MatPtAPMultEqual(A,B,PtAP,10,&flg);CHKERRQ(ierr); 417 if (!flg) { 418 ierr = PetscPrintf(PETSC_COMM_WORLD,"Error with PtAP\n");CHKERRQ(ierr); 419 ierr = MatMatMult(A,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&T);CHKERRQ(ierr); 420 ierr = MatTransposeMatMult(B,T,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&T2);CHKERRQ(ierr); 421 ierr = MatAXPY(T2,-1.0,PtAP,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 422 ierr = MatView(T2,NULL);CHKERRQ(ierr); 423 ierr = MatDestroy(&T2);CHKERRQ(ierr); 424 ierr = MatDestroy(&T);CHKERRQ(ierr); 425 } 426 ierr = PetscMalloc1((k+ldr)*M,&dataR);CHKERRQ(ierr); 427 ierr = MatCreateDense(PETSC_COMM_WORLD,PETSC_DECIDE,m,K,M,dataR,&R);CHKERRQ(ierr); 428 ierr = MatDenseSetLDA(R,k+ldr);CHKERRQ(ierr); 429 ierr = MatSetRandom(R,NULL);CHKERRQ(ierr); 430 if (testrart) { /* fails for AIJCUSPARSE because RA operation is not defined */ 431 ierr = MatRARt(A,R,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&RARt);CHKERRQ(ierr); 432 ierr = MatRARtMultEqual(A,R,RARt,10,&flg);CHKERRQ(ierr); 433 if (!flg) { 434 ierr = PetscPrintf(PETSC_COMM_WORLD,"Error with RARt\n");CHKERRQ(ierr); 435 ierr = MatMatTransposeMult(A,R,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&T);CHKERRQ(ierr); 436 ierr = MatMatMult(R,T,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&T2);CHKERRQ(ierr); 437 ierr = MatAXPY(T2,-1.0,RARt,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 438 ierr = MatView(T2,NULL);CHKERRQ(ierr); 439 ierr = MatDestroy(&T2);CHKERRQ(ierr); 440 ierr = MatDestroy(&T);CHKERRQ(ierr); 441 } 442 } 443 } 444 445 /* Test MatDenseGetColumnVec and friends */ 446 ierr = MatMatMult(A,B,MAT_REUSE_MATRIX,PETSC_DEFAULT,&X);CHKERRQ(ierr); 447 ierr = MatMatMult(A,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&T);CHKERRQ(ierr); 448 ierr = MatDuplicate(T,MAT_DO_NOT_COPY_VALUES,&T2);CHKERRQ(ierr); 449 for (k=0;k<K;k++) { 450 Vec Xv,Tv,T2v; 451 452 ierr = MatDenseGetColumnVecRead(X,k,&Xv);CHKERRQ(ierr); 453 ierr = MatDenseGetColumnVec(T,k,&Tv);CHKERRQ(ierr); 454 ierr = MatDenseGetColumnVecWrite(T2,k,&T2v);CHKERRQ(ierr); 455 ierr = VecCopy(Xv,T2v);CHKERRQ(ierr); 456 ierr = VecAXPY(Tv,-1.,Xv);CHKERRQ(ierr); 457 ierr = MatDenseRestoreColumnVecRead(X,k,&Xv);CHKERRQ(ierr); 458 ierr = MatDenseRestoreColumnVec(T,k,&Tv);CHKERRQ(ierr); 459 ierr = MatDenseRestoreColumnVecWrite(T2,k,&T2v);CHKERRQ(ierr); 460 } 461 ierr = MatNorm(T,NORM_FROBENIUS,&err);CHKERRQ(ierr); 462 if (err > PETSC_SMALL) { 463 ierr = PetscPrintf(PETSC_COMM_WORLD,"Error with MatDenseGetColumnVec\n");CHKERRQ(ierr); 464 ierr = MatView(T,NULL);CHKERRQ(ierr); 465 } 466 ierr = MatAXPY(T2,-1.,X,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 467 ierr = MatNorm(T2,NORM_FROBENIUS,&err);CHKERRQ(ierr); 468 if (err > PETSC_SMALL) { 469 ierr = PetscPrintf(PETSC_COMM_WORLD,"Error with MatDenseGetColumnVecWrite\n");CHKERRQ(ierr); 470 ierr = MatView(T2,NULL);CHKERRQ(ierr); 471 } 472 ierr = MatDestroy(&T);CHKERRQ(ierr); 473 ierr = MatDestroy(&T2);CHKERRQ(ierr); 474 475 /* Test with MatShell */ 476 ierr = MatDuplicate(A,MAT_COPY_VALUES,&T);CHKERRQ(ierr); 477 ierr = MatConvert(T,MATSHELL,MAT_INITIAL_MATRIX,&T2);CHKERRQ(ierr); 478 ierr = MatDestroy(&T);CHKERRQ(ierr); 479 480 /* scale matrix */ 481 ierr = MatScale(A,2.0);CHKERRQ(ierr); 482 ierr = MatScale(T2,2.0);CHKERRQ(ierr); 483 ierr = MatCreateVecs(A,&r,&l);CHKERRQ(ierr); 484 ierr = VecSetRandom(r,NULL);CHKERRQ(ierr); 485 ierr = VecSetRandom(l,NULL);CHKERRQ(ierr); 486 ierr = MatCreateVecs(T2,&rs,&ls);CHKERRQ(ierr); 487 ierr = VecCopy(r,rs);CHKERRQ(ierr); 488 ierr = VecCopy(l,ls);CHKERRQ(ierr); 489 if (testproj) { 490 ierr = MatDiagonalScale(A,r,r);CHKERRQ(ierr); 491 ierr = MatDiagonalScale(T2,rs,rs);CHKERRQ(ierr); 492 } else { 493 ierr = MatDiagonalScale(A,l,r);CHKERRQ(ierr); 494 ierr = MatDiagonalScale(T2,ls,rs);CHKERRQ(ierr); 495 } 496 ierr = MatDuplicate(A,MAT_COPY_VALUES,&T);CHKERRQ(ierr); 497 ierr = MatAXPY(A,4.5,T,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 498 ierr = MatAXPY(T2,4.5,T,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr); 499 ierr = MatMultEqual(T2,A,10,&flg);CHKERRQ(ierr); 500 if (!flg) { 501 ierr = PetscPrintf(PETSC_COMM_WORLD,"Error with MATSHELL (MatMult)\n");CHKERRQ(ierr); 502 } 503 ierr = MatMultTransposeEqual(T2,A,10,&flg);CHKERRQ(ierr); 504 if (!flg) { 505 ierr = PetscPrintf(PETSC_COMM_WORLD,"Error with MATSHELL (MatMultTranspose)\n");CHKERRQ(ierr); 506 } 507 ierr = MatDestroy(&T);CHKERRQ(ierr); 508 ierr = VecDestroy(&ls);CHKERRQ(ierr); 509 ierr = VecDestroy(&rs);CHKERRQ(ierr); 510 ierr = VecDestroy(&l);CHKERRQ(ierr); 511 ierr = VecDestroy(&r);CHKERRQ(ierr); 512 513 /* recompute projections, test reusage */ 514 if (PtAP) { ierr = MatPtAP(A,B,MAT_REUSE_MATRIX,PETSC_DEFAULT,&PtAP);CHKERRQ(ierr); } 515 if (RARt) { ierr = MatRARt(A,R,MAT_REUSE_MATRIX,PETSC_DEFAULT,&RARt);CHKERRQ(ierr); } 516 if (testshellops) { /* test callbacks for user defined MatProducts */ 517 ierr = MatShellSetMatProductOperation(T2,MATPRODUCT_AB,NULL,MyMatShellMatMultNumeric,NULL,MATDENSE,MATDENSE);CHKERRQ(ierr); 518 ierr = MatShellSetMatProductOperation(T2,MATPRODUCT_AB,NULL,MyMatShellMatMultNumeric,NULL,MATDENSECUDA,MATDENSECUDA);CHKERRQ(ierr); 519 ierr = MatShellSetMatProductOperation(T2,MATPRODUCT_AtB,NULL,MyMatTransposeShellMatMultNumeric,NULL,MATDENSE,MATDENSE);CHKERRQ(ierr); 520 ierr = MatShellSetMatProductOperation(T2,MATPRODUCT_AtB,NULL,MyMatTransposeShellMatMultNumeric,NULL,MATDENSECUDA,MATDENSECUDA);CHKERRQ(ierr); 521 ierr = MatShellSetMatProductOperation(T2,MATPRODUCT_ABt,NULL,MyMatShellMatTransposeMultNumeric,NULL,MATDENSE,MATDENSE);CHKERRQ(ierr); 522 ierr = MatShellSetMatProductOperation(T2,MATPRODUCT_ABt,NULL,MyMatShellMatTransposeMultNumeric,NULL,MATDENSECUDA,MATDENSECUDA);CHKERRQ(ierr); 523 if (testproj) { 524 ierr = MatShellSetMatProductOperation(T2,MATPRODUCT_PtAP,MyPtShellPMultSymbolic,MyPtShellPMultNumeric,proj_destroy,MATDENSE,MATSHELL);CHKERRQ(ierr); 525 ierr = MatShellSetMatProductOperation(T2,MATPRODUCT_PtAP,MyPtShellPMultSymbolic,MyPtShellPMultNumeric,proj_destroy,MATDENSECUDA,MATSHELL);CHKERRQ(ierr); 526 ierr = MatShellSetMatProductOperation(T2,MATPRODUCT_RARt,MyRShellRtMultSymbolic,MyRShellRtMultNumeric,proj_destroy,MATDENSE,MATSHELL);CHKERRQ(ierr); 527 ierr = MatShellSetMatProductOperation(T2,MATPRODUCT_RARt,MyRShellRtMultSymbolic,MyRShellRtMultNumeric,proj_destroy,MATDENSECUDA,MATSHELL);CHKERRQ(ierr); 528 } 529 } 530 ierr = CheckLocal(B,X,aB,aX);CHKERRQ(ierr); 531 /* we either use the shell operations or the loop over columns code, applying the operator */ 532 ierr = MatMatMult(T2,B,MAT_REUSE_MATRIX,PETSC_DEFAULT,&X);CHKERRQ(ierr); 533 ierr = CheckLocal(B,X,aB,aX);CHKERRQ(ierr); 534 ierr = MatMatMultEqual(T2,B,X,10,&flg);CHKERRQ(ierr); 535 if (!flg) { 536 ierr = PetscPrintf(PETSC_COMM_WORLD,"Error with reusage (MATSHELL)\n");CHKERRQ(ierr); 537 ierr = MatMatMult(A,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&T);CHKERRQ(ierr); 538 ierr = MatAXPY(T,-1.0,X,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 539 ierr = MatView(T,NULL);CHKERRQ(ierr); 540 ierr = MatDestroy(&T);CHKERRQ(ierr); 541 } 542 if (testproj) { 543 ierr = MatPtAPMultEqual(T2,B,PtAP,10,&flg);CHKERRQ(ierr); 544 if (!flg) { 545 ierr = PetscPrintf(PETSC_COMM_WORLD,"Error with PtAP\n");CHKERRQ(ierr); 546 } 547 if (testshellops) { /* projections fail if the product operations are not specified */ 548 ierr = MatPtAP(T2,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&T);CHKERRQ(ierr); 549 ierr = MatPtAP(T2,B,MAT_REUSE_MATRIX,PETSC_DEFAULT,&T);CHKERRQ(ierr); 550 ierr = MatPtAPMultEqual(T2,B,T,10,&flg);CHKERRQ(ierr); 551 if (!flg) { 552 Mat TE; 553 554 ierr = PetscPrintf(PETSC_COMM_WORLD,"Error with PtAP (user defined)\n");CHKERRQ(ierr); 555 ierr = MatComputeOperator(T,MATDENSE,&TE);CHKERRQ(ierr); 556 ierr = MatView(TE,NULL);CHKERRQ(ierr); 557 ierr = MatView(PtAP,NULL);CHKERRQ(ierr); 558 ierr = MatAXPY(TE,-1.0,PtAP,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 559 ierr = MatView(TE,NULL);CHKERRQ(ierr); 560 ierr = MatDestroy(&TE);CHKERRQ(ierr); 561 } 562 ierr = MatDestroy(&T);CHKERRQ(ierr); 563 } 564 if (RARt) { 565 ierr = MatRARtMultEqual(T2,R,RARt,10,&flg);CHKERRQ(ierr); 566 if (!flg) { 567 ierr = PetscPrintf(PETSC_COMM_WORLD,"Error with RARt\n");CHKERRQ(ierr); 568 } 569 } 570 if (testshellops) { 571 ierr = MatRARt(T2,R,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&T);CHKERRQ(ierr); 572 ierr = MatRARt(T2,R,MAT_REUSE_MATRIX,PETSC_DEFAULT,&T);CHKERRQ(ierr); 573 ierr = MatRARtMultEqual(T2,R,T,10,&flg);CHKERRQ(ierr); 574 if (!flg) { 575 Mat TE; 576 577 ierr = PetscPrintf(PETSC_COMM_WORLD,"Error with RARt (user defined)\n");CHKERRQ(ierr); 578 ierr = MatComputeOperator(T,MATDENSE,&TE);CHKERRQ(ierr); 579 ierr = MatView(TE,NULL);CHKERRQ(ierr); 580 if (RARt) { 581 ierr = MatView(RARt,NULL);CHKERRQ(ierr); 582 ierr = MatAXPY(TE,-1.0,RARt,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 583 ierr = MatView(TE,NULL);CHKERRQ(ierr); 584 } 585 ierr = MatDestroy(&TE);CHKERRQ(ierr); 586 } 587 ierr = MatDestroy(&T);CHKERRQ(ierr); 588 } 589 } 590 591 if (testmattmat) { /* we either use the shell operations or the loop over columns code applying the transposed operator */ 592 ierr = MatTransposeMatMult(T2,X,MAT_REUSE_MATRIX,PETSC_DEFAULT,&B);CHKERRQ(ierr); 593 ierr = CheckLocal(B,X,aB,aX);CHKERRQ(ierr); 594 ierr = MatTransposeMatMultEqual(T2,X,B,10,&flg);CHKERRQ(ierr); 595 if (!flg) { 596 ierr = PetscPrintf(PETSC_COMM_WORLD,"Error with reusage (MatTranspose, MATSHELL)\n");CHKERRQ(ierr); 597 ierr = MatTransposeMatMult(A,X,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&T);CHKERRQ(ierr); 598 ierr = MatAXPY(T,-1.0,B,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 599 ierr = MatView(T,NULL);CHKERRQ(ierr); 600 ierr = MatDestroy(&T);CHKERRQ(ierr); 601 } 602 } 603 if (testmatmatt && testshellops) { /* only when shell operations are set */ 604 ierr = MatMatTransposeMult(T2,Bt,MAT_REUSE_MATRIX,PETSC_DEFAULT,&X);CHKERRQ(ierr); 605 ierr = CheckLocal(Bt,X,aBt,aX);CHKERRQ(ierr); 606 ierr = MatMatTransposeMultEqual(T2,Bt,X,10,&flg);CHKERRQ(ierr); 607 if (!flg) { 608 ierr = PetscPrintf(PETSC_COMM_WORLD,"Error with reusage (MatMatTranspose, MATSHELL)\n");CHKERRQ(ierr); 609 ierr = MatMatTransposeMult(A,Bt,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&T);CHKERRQ(ierr); 610 ierr = MatAXPY(T,-1.0,X,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 611 ierr = MatView(T,NULL);CHKERRQ(ierr); 612 ierr = MatDestroy(&T);CHKERRQ(ierr); 613 } 614 } 615 ierr = MatDestroy(&T2);CHKERRQ(ierr); 616 617 if (testnest) { /* test with MatNest */ 618 Mat NA; 619 const char *vtype; 620 621 ierr = MatCreateNest(PETSC_COMM_WORLD,1,NULL,1,NULL,&A,&NA);CHKERRQ(ierr); 622 /* needed to test against CUSPARSE matrices */ 623 ierr = MatGetVecType(A,&vtype);CHKERRQ(ierr); 624 ierr = MatSetVecType(NA,vtype);CHKERRQ(ierr); 625 ierr = MatViewFromOptions(NA,NULL,"-NA_view");CHKERRQ(ierr); 626 ierr = MatMatMult(NA,B,MAT_REUSE_MATRIX,PETSC_DEFAULT,&X);CHKERRQ(ierr); 627 ierr = CheckLocal(B,X,aB,aX);CHKERRQ(ierr); 628 ierr = MatMatMultEqual(NA,B,X,10,&flg);CHKERRQ(ierr); 629 if (!flg) { 630 ierr = PetscPrintf(PETSC_COMM_WORLD,"Error with Nest\n");CHKERRQ(ierr); 631 ierr = MatMatMult(NA,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&T);CHKERRQ(ierr); 632 ierr = MatAXPY(T,-1.0,X,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 633 ierr = MatView(T,NULL);CHKERRQ(ierr); 634 ierr = MatDestroy(&T);CHKERRQ(ierr); 635 } 636 ierr = MatDestroy(&NA);CHKERRQ(ierr); 637 } 638 639 if (testtranspose) { /* test with Transpose */ 640 Mat TA; 641 642 ierr = MatCreateTranspose(A,&TA);CHKERRQ(ierr); 643 ierr = MatMatMult(TA,X,MAT_REUSE_MATRIX,PETSC_DEFAULT,&B);CHKERRQ(ierr); 644 ierr = CheckLocal(B,X,aB,aX);CHKERRQ(ierr); 645 ierr = MatMatMultEqual(TA,X,B,10,&flg);CHKERRQ(ierr); 646 if (!flg) { 647 ierr = PetscPrintf(PETSC_COMM_WORLD,"Error with Transpose\n");CHKERRQ(ierr); 648 ierr = MatMatMult(TA,X,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&T);CHKERRQ(ierr); 649 ierr = MatAXPY(T,-1.0,B,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 650 ierr = MatView(T,NULL);CHKERRQ(ierr); 651 ierr = MatDestroy(&T);CHKERRQ(ierr); 652 } 653 ierr = MatDestroy(&TA);CHKERRQ(ierr); 654 } 655 656 if (testhtranspose) { /* test with Hermitian Transpose */ 657 Mat TA; 658 659 ierr = MatCreateHermitianTranspose(A,&TA);CHKERRQ(ierr); 660 ierr = MatMatMult(TA,X,MAT_REUSE_MATRIX,PETSC_DEFAULT,&B);CHKERRQ(ierr); 661 ierr = CheckLocal(B,X,aB,aX);CHKERRQ(ierr); 662 ierr = MatMatMultEqual(TA,X,B,10,&flg);CHKERRQ(ierr); 663 if (!flg) { 664 ierr = PetscPrintf(PETSC_COMM_WORLD,"Error with Transpose\n");CHKERRQ(ierr); 665 ierr = MatMatMult(TA,X,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&T);CHKERRQ(ierr); 666 ierr = MatAXPY(T,-1.0,B,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 667 ierr = MatView(T,NULL);CHKERRQ(ierr); 668 ierr = MatDestroy(&T);CHKERRQ(ierr); 669 } 670 ierr = MatDestroy(&TA);CHKERRQ(ierr); 671 } 672 673 if (testtt) { /* test with Transpose(Transpose) */ 674 Mat TA, TTA; 675 676 ierr = MatCreateTranspose(A,&TA);CHKERRQ(ierr); 677 ierr = MatCreateTranspose(TA,&TTA);CHKERRQ(ierr); 678 ierr = MatMatMult(TTA,B,MAT_REUSE_MATRIX,PETSC_DEFAULT,&X);CHKERRQ(ierr); 679 ierr = CheckLocal(B,X,aB,aX);CHKERRQ(ierr); 680 ierr = MatMatMultEqual(TTA,B,X,10,&flg);CHKERRQ(ierr); 681 if (!flg) { 682 ierr = PetscPrintf(PETSC_COMM_WORLD,"Error with Transpose(Transpose)\n");CHKERRQ(ierr); 683 ierr = MatMatMult(TTA,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&T);CHKERRQ(ierr); 684 ierr = MatAXPY(T,-1.0,X,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 685 ierr = MatView(T,NULL);CHKERRQ(ierr); 686 ierr = MatDestroy(&T);CHKERRQ(ierr); 687 } 688 ierr = MatDestroy(&TA);CHKERRQ(ierr); 689 ierr = MatDestroy(&TTA);CHKERRQ(ierr); 690 } 691 692 if (testcircular) { /* test circular */ 693 Mat AB; 694 695 ierr = MatMatMult(A,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&AB);CHKERRQ(ierr); 696 ierr = MatMatMult(A,B,MAT_REUSE_MATRIX,PETSC_DEFAULT,&X);CHKERRQ(ierr); 697 ierr = CheckLocal(B,X,aB,aX);CHKERRQ(ierr); 698 if (M == N && N == K) { 699 ierr = MatMatMult(A,X,MAT_REUSE_MATRIX,PETSC_DEFAULT,&B);CHKERRQ(ierr); 700 } else { 701 ierr = MatTransposeMatMult(A,X,MAT_REUSE_MATRIX,PETSC_DEFAULT,&B);CHKERRQ(ierr); 702 } 703 ierr = CheckLocal(B,X,aB,aX);CHKERRQ(ierr); 704 ierr = MatDestroy(&AB);CHKERRQ(ierr); 705 } 706 ierr = MatDestroy(&X);CHKERRQ(ierr); 707 ierr = MatDestroy(&Bt);CHKERRQ(ierr); 708 ierr = MatDestroy(&B);CHKERRQ(ierr); 709 ierr = MatDestroy(&A);CHKERRQ(ierr); 710 ierr = MatDestroy(&R);CHKERRQ(ierr); 711 ierr = MatDestroy(&PtAP);CHKERRQ(ierr); 712 ierr = MatDestroy(&RARt);CHKERRQ(ierr); 713 ierr = PetscFree(dataX);CHKERRQ(ierr); 714 ierr = PetscFree(dataB);CHKERRQ(ierr); 715 ierr = PetscFree(dataR);CHKERRQ(ierr); 716 ierr = PetscFree(dataBt);CHKERRQ(ierr); 717 ierr = PetscFinalize(); 718 return ierr; 719 } 720 721 /*TEST 722 723 test: 724 suffix: 1 725 args: -local {{0 1}} -testshellops 726 727 test: 728 output_file: output/ex70_1.out 729 requires: cuda 730 suffix: 1_cuda 731 args: -local {{0 1}} -xgpu {{0 1}} -bgpu {{0 1}} -A_mat_type {{seqaijcusparse seqaij}} -testshellops {{0 1}} 732 733 test: 734 output_file: output/ex70_1.out 735 nsize: 2 736 suffix: 1_par 737 args: -local {{0 1}} -testmatmatt 0 738 739 test: 740 output_file: output/ex70_1.out 741 requires: cuda 742 nsize: 2 743 suffix: 1_par_cuda 744 args: -local {{0 1}} -xgpu {{0 1}} -bgpu {{0 1}} -A_mat_type {{mpiaijcusparse mpiaij}} -testnest 0 -testmatmatt 0 -matmatmult_Bbn 3 745 746 test: 747 output_file: output/ex70_1.out 748 suffix: 2 749 nsize: 1 750 args: -M {{7 11}} -N {{12 9}} -K {{1 3}} -local {{0 1}} 751 752 testset: 753 requires: cuda 754 output_file: output/ex70_1.out 755 nsize: 1 756 args: -M 7 -N 9 -K 2 -local {{0 1}} -testnest 0 -A_mat_type {{seqdensecuda seqdense}} -xgpu {{0 1}} -bgpu {{0 1}} 757 test: 758 requires: !complex 759 suffix: 2_cuda_real 760 test: 761 # complex+single gives a little bigger error in the MatDenseGetColumnVec test 762 requires: complex !single 763 suffix: 2_cuda_complex 764 765 test: 766 output_file: output/ex70_1.out 767 suffix: 2_par 768 nsize: 2 769 args: -M {{7 11}} -N {{12 9}} -K {{1 3}} -local {{0 1}} -testcircular -testmatmatt 0 770 771 test: 772 requires: cuda 773 output_file: output/ex70_1.out 774 suffix: 2_par_cuda 775 nsize: 2 776 args: -M 11 -N 9 -K 1 -local {{0 1}} -testcircular 0 -A_mat_type mpiaijcusparse -xgpu -bgpu -testnest 0 -testmatmatt 0 777 778 test: 779 output_file: output/ex70_1.out 780 suffix: 3 781 nsize: {{1 3}} 782 args: -M 13 -N 13 -K {{1 3}} -local {{0 1}} -A_mat_type sbaij -symm -testproj 0 -testmatmatt 0 783 784 test: 785 output_file: output/ex70_1.out 786 suffix: 4 787 nsize: 1 788 args: -M 3 -N 3 -K 3 -local {{0 1}} -testcircular 789 790 test: 791 output_file: output/ex70_1.out 792 suffix: 5 793 nsize: {{2 4}} 794 args: -M 3 -N 3 -K 3 -local {{0 1}} -testcircular -testmatmatt 0 795 796 test: 797 output_file: output/ex70_1.out 798 suffix: 6 799 nsize: 1 800 args: -M {{1 3}} -N {{2 5}} -K {{1 2}} -local {{0 1}} -testcircular 801 802 test: 803 output_file: output/ex70_1.out 804 suffix: 7 805 nsize: 1 806 args: -M 13 -N 13 -K {{1 3}} -local {{0 1}} -A_mat_type dense -testnest -testcircular 807 808 TEST*/ 809