1 2 static char help[] = "Tests sequential and parallel MatMatMult() and MatPtAP(), MatTransposeMatMult(), sequential MatMatTransposeMult(), MatRARt()\n\ 3 Input arguments are:\n\ 4 -f0 <input_file> -f1 <input_file> -f2 <input_file> -f3 <input_file> : file to load\n\n"; 5 /* Example of usage: 6 ./ex94 -f0 <A_binary> -f1 <B_binary> -matmatmult_mat_view ascii::ascii_info -matmatmulttr_mat_view 7 mpiexec -n 3 ./ex94 -f0 medium -f1 medium -f2 arco1 -f3 arco1 -matmatmult_mat_view 8 */ 9 10 #include <petscmat.h> 11 12 /* 13 B = A - B 14 norm = norm(B) 15 */ 16 PetscErrorCode MatNormDifference(Mat A,Mat B,PetscReal *norm) 17 { 18 PetscFunctionBegin; 19 PetscCall(MatAXPY(B,-1.0,A,DIFFERENT_NONZERO_PATTERN)); 20 PetscCall(MatNorm(B,NORM_FROBENIUS,norm)); 21 PetscFunctionReturn(0); 22 } 23 24 int main(int argc,char **args) 25 { 26 Mat A,A_save,B,AT,ATT,BT,BTT,P,R,C,C1; 27 Vec x,v1,v2,v3,v4; 28 PetscViewer viewer; 29 PetscMPIInt size,rank; 30 PetscInt i,m,n,j,*idxn,M,N,nzp,rstart,rend; 31 PetscReal norm,norm_abs,norm_tmp,fill=4.0; 32 PetscRandom rdm; 33 char file[4][128]; 34 PetscBool flg,preload = PETSC_TRUE; 35 PetscScalar *a,rval,alpha,none = -1.0; 36 PetscBool Test_MatMatMult=PETSC_TRUE,Test_MatMatTr=PETSC_TRUE,Test_MatPtAP=PETSC_TRUE,Test_MatRARt=PETSC_TRUE,Test_MatMatMatMult=PETSC_TRUE; 37 PetscBool Test_MatAXPY=PETSC_FALSE,view=PETSC_FALSE; 38 PetscInt pm,pn,pM,pN; 39 MatInfo info; 40 PetscBool seqaij; 41 MatType mattype; 42 Mat Cdensetest,Pdense,Cdense,Adense; 43 44 PetscFunctionBeginUser; 45 PetscCall(PetscInitialize(&argc,&args,(char*)0,help)); 46 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 47 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank)); 48 49 PetscCall(PetscOptionsGetReal(NULL,NULL,"-fill",&fill,NULL)); 50 PetscCall(PetscOptionsGetBool(NULL,NULL,"-matops_view",&view,NULL)); 51 if (view) PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO)); 52 53 /* Load the matrices A_save and B */ 54 PetscCall(PetscOptionsGetString(NULL,NULL,"-f0",file[0],sizeof(file[0]),&flg)); 55 PetscCheck(flg,PETSC_COMM_WORLD,PETSC_ERR_USER,"Must indicate a file name for small matrix A with the -f0 option."); 56 PetscCall(PetscOptionsGetString(NULL,NULL,"-f1",file[1],sizeof(file[1]),&flg)); 57 PetscCheck(flg,PETSC_COMM_WORLD,PETSC_ERR_USER,"Must indicate a file name for small matrix B with the -f1 option."); 58 PetscCall(PetscOptionsGetString(NULL,NULL,"-f2",file[2],sizeof(file[2]),&flg)); 59 if (!flg) { 60 preload = PETSC_FALSE; 61 } else { 62 PetscCall(PetscOptionsGetString(NULL,NULL,"-f3",file[3],sizeof(file[3]),&flg)); 63 PetscCheck(flg,PETSC_COMM_WORLD,PETSC_ERR_USER,"Must indicate a file name for test matrix B with the -f3 option."); 64 } 65 66 PetscPreLoadBegin(preload,"Load system"); 67 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[2*PetscPreLoadIt],FILE_MODE_READ,&viewer)); 68 PetscCall(MatCreate(PETSC_COMM_WORLD,&A_save)); 69 PetscCall(MatSetFromOptions(A_save)); 70 PetscCall(MatLoad(A_save,viewer)); 71 PetscCall(PetscViewerDestroy(&viewer)); 72 73 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[2*PetscPreLoadIt+1],FILE_MODE_READ,&viewer)); 74 PetscCall(MatCreate(PETSC_COMM_WORLD,&B)); 75 PetscCall(MatSetFromOptions(B)); 76 PetscCall(MatLoad(B,viewer)); 77 PetscCall(PetscViewerDestroy(&viewer)); 78 79 PetscCall(MatGetType(B,&mattype)); 80 81 PetscCall(MatGetSize(B,&M,&N)); 82 nzp = PetscMax((PetscInt)(0.1*M),5); 83 PetscCall(PetscMalloc((nzp+1)*(sizeof(PetscInt)+sizeof(PetscScalar)),&idxn)); 84 a = (PetscScalar*)(idxn + nzp); 85 86 /* Create vectors v1 and v2 that are compatible with A_save */ 87 PetscCall(VecCreate(PETSC_COMM_WORLD,&v1)); 88 PetscCall(MatGetLocalSize(A_save,&m,NULL)); 89 PetscCall(VecSetSizes(v1,m,PETSC_DECIDE)); 90 PetscCall(VecSetFromOptions(v1)); 91 PetscCall(VecDuplicate(v1,&v2)); 92 93 PetscCall(PetscRandomCreate(PETSC_COMM_WORLD,&rdm)); 94 PetscCall(PetscRandomSetFromOptions(rdm)); 95 PetscCall(PetscOptionsGetReal(NULL,NULL,"-fill",&fill,NULL)); 96 97 /* Test MatAXPY() */ 98 /*-------------------*/ 99 PetscCall(PetscOptionsHasName(NULL,NULL,"-test_MatAXPY",&Test_MatAXPY)); 100 if (Test_MatAXPY) { 101 Mat Btmp; 102 PetscCall(MatDuplicate(A_save,MAT_COPY_VALUES,&A)); 103 PetscCall(MatDuplicate(B,MAT_COPY_VALUES,&Btmp)); 104 PetscCall(MatAXPY(A,-1.0,B,DIFFERENT_NONZERO_PATTERN)); /* A = -B + A_save */ 105 106 PetscCall(MatScale(A,-1.0)); /* A = -A = B - A_save */ 107 PetscCall(MatAXPY(Btmp,-1.0,A,DIFFERENT_NONZERO_PATTERN)); /* Btmp = -A + B = A_save */ 108 PetscCall(MatMultEqual(A_save,Btmp,10,&flg)); 109 PetscCheck(flg,PETSC_COMM_SELF,PETSC_ERR_PLIB,"MatAXPY() is incorrect"); 110 PetscCall(MatDestroy(&A)); 111 PetscCall(MatDestroy(&Btmp)); 112 113 Test_MatMatMult = PETSC_FALSE; 114 Test_MatMatTr = PETSC_FALSE; 115 Test_MatPtAP = PETSC_FALSE; 116 Test_MatRARt = PETSC_FALSE; 117 Test_MatMatMatMult = PETSC_FALSE; 118 } 119 120 /* 1) Test MatMatMult() */ 121 /* ---------------------*/ 122 if (Test_MatMatMult) { 123 PetscCall(MatDuplicate(A_save,MAT_COPY_VALUES,&A)); 124 PetscCall(MatCreateTranspose(A,&AT)); 125 PetscCall(MatCreateTranspose(AT,&ATT)); 126 PetscCall(MatCreateTranspose(B,&BT)); 127 PetscCall(MatCreateTranspose(BT,&BTT)); 128 129 PetscCall(MatMatMult(AT,B,MAT_INITIAL_MATRIX,fill,&C)); 130 PetscCall(MatMatMultEqual(AT,B,C,10,&flg)); 131 PetscCheck(flg,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error: MatMatMult() for C=AT*B"); 132 PetscCall(MatDestroy(&C)); 133 134 PetscCall(MatMatMult(ATT,B,MAT_INITIAL_MATRIX,fill,&C)); 135 PetscCall(MatMatMultEqual(ATT,B,C,10,&flg)); 136 PetscCheck(flg,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error: MatMatMult() for C=ATT*B"); 137 PetscCall(MatDestroy(&C)); 138 139 PetscCall(MatMatMult(A,B,MAT_INITIAL_MATRIX,fill,&C)); 140 PetscCall(MatMatMultEqual(A,B,C,10,&flg)); 141 PetscCheck(flg,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error: MatMatMult() for reuse C=A*B"); 142 /* ATT has different matrix type as A (although they have same internal data structure), 143 we cannot call MatProductReplaceMats(ATT,NULL,NULL,C) and MatMatMult(ATT,B,MAT_REUSE_MATRIX,fill,&C) */ 144 PetscCall(MatDestroy(&C)); 145 146 PetscCall(MatMatMult(A,BTT,MAT_INITIAL_MATRIX,fill,&C)); 147 PetscCall(MatMatMultEqual(A,BTT,C,10,&flg)); 148 PetscCheck(flg,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error: MatMatMult() for C=A*BTT"); 149 PetscCall(MatDestroy(&C)); 150 151 PetscCall(MatMatMult(ATT,BTT,MAT_INITIAL_MATRIX,fill,&C)); 152 PetscCall(MatMatMultEqual(A,B,C,10,&flg)); 153 PetscCheck(flg,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error: MatMatMult()"); 154 PetscCall(MatDestroy(&C)); 155 156 PetscCall(MatDestroy(&BTT)); 157 PetscCall(MatDestroy(&BT)); 158 PetscCall(MatDestroy(&ATT)); 159 PetscCall(MatDestroy(&AT)); 160 161 PetscCall(MatMatMult(A,B,MAT_INITIAL_MATRIX,fill,&C)); 162 PetscCall(MatSetOptionsPrefix(C,"matmatmult_")); /* enable option '-matmatmult_' for matrix C */ 163 PetscCall(MatGetInfo(C,MAT_GLOBAL_SUM,&info)); 164 165 /* Test MAT_REUSE_MATRIX - reuse symbolic C */ 166 alpha=1.0; 167 for (i=0; i<2; i++) { 168 alpha -=0.1; 169 PetscCall(MatScale(A,alpha)); 170 PetscCall(MatMatMult(A,B,MAT_REUSE_MATRIX,fill,&C)); 171 } 172 PetscCall(MatMatMultEqual(A,B,C,10,&flg)); 173 PetscCheck(flg,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error: MatMatMult()"); 174 PetscCall(MatDestroy(&A)); 175 176 /* Test MatDuplicate() of C=A*B */ 177 PetscCall(MatDuplicate(C,MAT_COPY_VALUES,&C1)); 178 PetscCall(MatDestroy(&C1)); 179 PetscCall(MatDestroy(&C)); 180 } /* if (Test_MatMatMult) */ 181 182 /* 2) Test MatTransposeMatMult() and MatMatTransposeMult() */ 183 /* ------------------------------------------------------- */ 184 if (Test_MatMatTr) { 185 /* Create P */ 186 PetscInt PN,rstart,rend; 187 PN = M/2; 188 nzp = 5; /* num of nonzeros in each row of P */ 189 PetscCall(MatCreate(PETSC_COMM_WORLD,&P)); 190 PetscCall(MatSetSizes(P,PETSC_DECIDE,PETSC_DECIDE,M,PN)); 191 PetscCall(MatSetType(P,mattype)); 192 PetscCall(MatSeqAIJSetPreallocation(P,nzp,NULL)); 193 PetscCall(MatMPIAIJSetPreallocation(P,nzp,NULL,nzp,NULL)); 194 PetscCall(MatGetOwnershipRange(P,&rstart,&rend)); 195 for (i=0; i<nzp; i++) { 196 PetscCall(PetscRandomGetValue(rdm,&a[i])); 197 } 198 for (i=rstart; i<rend; i++) { 199 for (j=0; j<nzp; j++) { 200 PetscCall(PetscRandomGetValue(rdm,&rval)); 201 idxn[j] = (PetscInt)(PetscRealPart(rval)*PN); 202 } 203 PetscCall(MatSetValues(P,1,&i,nzp,idxn,a,ADD_VALUES)); 204 } 205 PetscCall(MatAssemblyBegin(P,MAT_FINAL_ASSEMBLY)); 206 PetscCall(MatAssemblyEnd(P,MAT_FINAL_ASSEMBLY)); 207 208 /* Create R = P^T */ 209 PetscCall(MatTranspose(P,MAT_INITIAL_MATRIX,&R)); 210 211 { /* Test R = P^T, C1 = R*B */ 212 PetscCall(MatMatMult(R,B,MAT_INITIAL_MATRIX,fill,&C1)); 213 PetscCall(MatTranspose(P,MAT_REUSE_MATRIX,&R)); 214 PetscCall(MatMatMult(R,B,MAT_REUSE_MATRIX,fill,&C1)); 215 PetscCall(MatDestroy(&C1)); 216 } 217 218 /* C = P^T*B */ 219 PetscCall(MatTransposeMatMult(P,B,MAT_INITIAL_MATRIX,fill,&C)); 220 PetscCall(MatGetInfo(C,MAT_GLOBAL_SUM,&info)); 221 222 /* Test MAT_REUSE_MATRIX - reuse symbolic C */ 223 PetscCall(MatTransposeMatMult(P,B,MAT_REUSE_MATRIX,fill,&C)); 224 if (view) { 225 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"C = P^T * B:\n")); 226 PetscCall(MatView(C,PETSC_VIEWER_STDOUT_WORLD)); 227 } 228 PetscCall(MatProductClear(C)); 229 if (view) { 230 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nC = P^T * B after MatProductClear():\n")); 231 PetscCall(MatView(C,PETSC_VIEWER_STDOUT_WORLD)); 232 } 233 234 /* Compare P^T*B and R*B */ 235 PetscCall(MatMatMult(R,B,MAT_INITIAL_MATRIX,fill,&C1)); 236 PetscCall(MatNormDifference(C,C1,&norm)); 237 PetscCheck(norm <= PETSC_SMALL,PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Error in MatTransposeMatMult(): %g",(double)norm); 238 PetscCall(MatDestroy(&C1)); 239 240 /* Test MatDuplicate() of C=P^T*B */ 241 PetscCall(MatDuplicate(C,MAT_COPY_VALUES,&C1)); 242 PetscCall(MatDestroy(&C1)); 243 PetscCall(MatDestroy(&C)); 244 245 /* C = B*R^T */ 246 PetscCall(PetscObjectTypeCompare((PetscObject)B,MATSEQAIJ,&seqaij)); 247 if (size == 1 && seqaij) { 248 PetscCall(MatMatTransposeMult(B,R,MAT_INITIAL_MATRIX,fill,&C)); 249 PetscCall(MatSetOptionsPrefix(C,"matmatmulttr_")); /* enable '-matmatmulttr_' for matrix C */ 250 PetscCall(MatGetInfo(C,MAT_GLOBAL_SUM,&info)); 251 252 /* Test MAT_REUSE_MATRIX - reuse symbolic C */ 253 PetscCall(MatMatTransposeMult(B,R,MAT_REUSE_MATRIX,fill,&C)); 254 255 /* Check */ 256 PetscCall(MatMatMult(B,P,MAT_INITIAL_MATRIX,fill,&C1)); 257 PetscCall(MatNormDifference(C,C1,&norm)); 258 PetscCheck(norm <= PETSC_SMALL,PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Error in MatMatTransposeMult() %g",(double)norm); 259 PetscCall(MatDestroy(&C1)); 260 PetscCall(MatDestroy(&C)); 261 } 262 PetscCall(MatDestroy(&P)); 263 PetscCall(MatDestroy(&R)); 264 } 265 266 /* 3) Test MatPtAP() */ 267 /*-------------------*/ 268 if (Test_MatPtAP) { 269 PetscInt PN; 270 Mat Cdup; 271 272 PetscCall(MatDuplicate(A_save,MAT_COPY_VALUES,&A)); 273 PetscCall(MatGetSize(A,&M,&N)); 274 PetscCall(MatGetLocalSize(A,&m,&n)); 275 276 PN = M/2; 277 nzp = (PetscInt)(0.1*PN+1); /* num of nozeros in each row of P */ 278 PetscCall(MatCreate(PETSC_COMM_WORLD,&P)); 279 PetscCall(MatSetSizes(P,PETSC_DECIDE,PETSC_DECIDE,N,PN)); 280 PetscCall(MatSetType(P,mattype)); 281 PetscCall(MatSeqAIJSetPreallocation(P,nzp,NULL)); 282 PetscCall(MatMPIAIJSetPreallocation(P,nzp,NULL,nzp,NULL)); 283 for (i=0; i<nzp; i++) { 284 PetscCall(PetscRandomGetValue(rdm,&a[i])); 285 } 286 PetscCall(MatGetOwnershipRange(P,&rstart,&rend)); 287 for (i=rstart; i<rend; i++) { 288 for (j=0; j<nzp; j++) { 289 PetscCall(PetscRandomGetValue(rdm,&rval)); 290 idxn[j] = (PetscInt)(PetscRealPart(rval)*PN); 291 } 292 PetscCall(MatSetValues(P,1,&i,nzp,idxn,a,ADD_VALUES)); 293 } 294 PetscCall(MatAssemblyBegin(P,MAT_FINAL_ASSEMBLY)); 295 PetscCall(MatAssemblyEnd(P,MAT_FINAL_ASSEMBLY)); 296 297 /* PetscCall(MatView(P,PETSC_VIEWER_STDOUT_WORLD)); */ 298 PetscCall(MatGetSize(P,&pM,&pN)); 299 PetscCall(MatGetLocalSize(P,&pm,&pn)); 300 PetscCall(MatPtAP(A,P,MAT_INITIAL_MATRIX,fill,&C)); 301 302 /* Test MAT_REUSE_MATRIX - reuse symbolic C */ 303 alpha=1.0; 304 for (i=0; i<2; i++) { 305 alpha -=0.1; 306 PetscCall(MatScale(A,alpha)); 307 PetscCall(MatPtAP(A,P,MAT_REUSE_MATRIX,fill,&C)); 308 } 309 310 /* Test PtAP ops with P Dense and A either AIJ or SeqDense (it assumes MatPtAP_XAIJ_XAIJ is fine) */ 311 PetscCall(PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&seqaij)); 312 if (seqaij) { 313 PetscCall(MatConvert(C,MATSEQDENSE,MAT_INITIAL_MATRIX,&Cdensetest)); 314 PetscCall(MatConvert(P,MATSEQDENSE,MAT_INITIAL_MATRIX,&Pdense)); 315 } else { 316 PetscCall(MatConvert(C,MATMPIDENSE,MAT_INITIAL_MATRIX,&Cdensetest)); 317 PetscCall(MatConvert(P,MATMPIDENSE,MAT_INITIAL_MATRIX,&Pdense)); 318 } 319 320 /* test with A(AIJ), Pdense -- call MatPtAP_Basic() when np>1 */ 321 PetscCall(MatPtAP(A,Pdense,MAT_INITIAL_MATRIX,fill,&Cdense)); 322 PetscCall(MatPtAP(A,Pdense,MAT_REUSE_MATRIX,fill,&Cdense)); 323 PetscCall(MatPtAPMultEqual(A,Pdense,Cdense,10,&flg)); 324 PetscCheck(flg,PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Error in MatPtAP with A AIJ and P Dense"); 325 PetscCall(MatDestroy(&Cdense)); 326 327 /* test with A SeqDense */ 328 if (seqaij) { 329 PetscCall(MatConvert(A,MATSEQDENSE,MAT_INITIAL_MATRIX,&Adense)); 330 PetscCall(MatPtAP(Adense,Pdense,MAT_INITIAL_MATRIX,fill,&Cdense)); 331 PetscCall(MatPtAP(Adense,Pdense,MAT_REUSE_MATRIX,fill,&Cdense)); 332 PetscCall(MatPtAPMultEqual(Adense,Pdense,Cdense,10,&flg)); 333 PetscCheck(flg,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error in MatPtAP with A SeqDense and P SeqDense"); 334 PetscCall(MatDestroy(&Cdense)); 335 PetscCall(MatDestroy(&Adense)); 336 } 337 PetscCall(MatDestroy(&Cdensetest)); 338 PetscCall(MatDestroy(&Pdense)); 339 340 /* Test MatDuplicate() of C=PtAP and MatView(Cdup,...) */ 341 PetscCall(MatDuplicate(C,MAT_COPY_VALUES,&Cdup)); 342 if (view) { 343 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nC = P^T * A * P:\n")); 344 PetscCall(MatView(C,PETSC_VIEWER_STDOUT_WORLD)); 345 346 PetscCall(MatProductClear(C)); 347 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nC = P^T * A * P after MatProductClear():\n")); 348 PetscCall(MatView(C,PETSC_VIEWER_STDOUT_WORLD)); 349 350 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nCdup:\n")); 351 PetscCall(MatView(Cdup,PETSC_VIEWER_STDOUT_WORLD)); 352 } 353 PetscCall(MatDestroy(&Cdup)); 354 355 if (size>1 || !seqaij) Test_MatRARt = PETSC_FALSE; 356 /* 4) Test MatRARt() */ 357 /* ----------------- */ 358 if (Test_MatRARt) { 359 Mat R, RARt, Rdense, RARtdense; 360 PetscCall(MatTranspose(P,MAT_INITIAL_MATRIX,&R)); 361 362 /* Test MatRARt_Basic(), MatMatMatMult_Basic() */ 363 PetscCall(MatConvert(R,MATDENSE,MAT_INITIAL_MATRIX,&Rdense)); 364 PetscCall(MatRARt(A,Rdense,MAT_INITIAL_MATRIX,2.0,&RARtdense)); 365 PetscCall(MatRARt(A,Rdense,MAT_REUSE_MATRIX,2.0,&RARtdense)); 366 367 PetscCall(MatConvert(RARtdense,MATAIJ,MAT_INITIAL_MATRIX,&RARt)); 368 PetscCall(MatNormDifference(C,RARt,&norm)); 369 PetscCheck(norm <= PETSC_SMALL,PETSC_COMM_SELF,PETSC_ERR_PLIB,"|PtAP - RARtdense| = %g",(double)norm); 370 PetscCall(MatDestroy(&Rdense)); 371 PetscCall(MatDestroy(&RARtdense)); 372 PetscCall(MatDestroy(&RARt)); 373 374 /* Test MatRARt() for aij matrices */ 375 PetscCall(MatRARt(A,R,MAT_INITIAL_MATRIX,2.0,&RARt)); 376 PetscCall(MatRARt(A,R,MAT_REUSE_MATRIX,2.0,&RARt)); 377 PetscCall(MatNormDifference(C,RARt,&norm)); 378 PetscCheck(norm <= PETSC_SMALL,PETSC_COMM_SELF,PETSC_ERR_PLIB,"|PtAP - RARt| = %g",(double)norm); 379 PetscCall(MatDestroy(&R)); 380 PetscCall(MatDestroy(&RARt)); 381 } 382 383 if (Test_MatMatMatMult && size == 1) { 384 Mat R, RAP; 385 PetscCall(MatTranspose(P,MAT_INITIAL_MATRIX,&R)); 386 PetscCall(MatMatMatMult(R,A,P,MAT_INITIAL_MATRIX,2.0,&RAP)); 387 PetscCall(MatMatMatMult(R,A,P,MAT_REUSE_MATRIX,2.0,&RAP)); 388 PetscCall(MatNormDifference(C,RAP,&norm)); 389 PetscCheck(norm <= PETSC_SMALL,PETSC_COMM_SELF,PETSC_ERR_PLIB,"PtAP != RAP %g",(double)norm); 390 PetscCall(MatDestroy(&R)); 391 PetscCall(MatDestroy(&RAP)); 392 } 393 394 /* Create vector x that is compatible with P */ 395 PetscCall(VecCreate(PETSC_COMM_WORLD,&x)); 396 PetscCall(MatGetLocalSize(P,&m,&n)); 397 PetscCall(VecSetSizes(x,n,PETSC_DECIDE)); 398 PetscCall(VecSetFromOptions(x)); 399 400 PetscCall(VecCreate(PETSC_COMM_WORLD,&v3)); 401 PetscCall(VecSetSizes(v3,n,PETSC_DECIDE)); 402 PetscCall(VecSetFromOptions(v3)); 403 PetscCall(VecDuplicate(v3,&v4)); 404 405 norm = 0.0; 406 for (i=0; i<10; i++) { 407 PetscCall(VecSetRandom(x,rdm)); 408 PetscCall(MatMult(P,x,v1)); 409 PetscCall(MatMult(A,v1,v2)); /* v2 = A*P*x */ 410 411 PetscCall(MatMultTranspose(P,v2,v3)); /* v3 = Pt*A*P*x */ 412 PetscCall(MatMult(C,x,v4)); /* v3 = C*x */ 413 PetscCall(VecNorm(v4,NORM_2,&norm_abs)); 414 PetscCall(VecAXPY(v4,none,v3)); 415 PetscCall(VecNorm(v4,NORM_2,&norm_tmp)); 416 417 norm_tmp /= norm_abs; 418 if (norm_tmp > norm) norm = norm_tmp; 419 } 420 PetscCheck(norm < PETSC_SMALL,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error: MatPtAP(), |v1 - v2|: %g",(double)norm); 421 422 PetscCall(MatDestroy(&A)); 423 PetscCall(MatDestroy(&P)); 424 PetscCall(MatDestroy(&C)); 425 PetscCall(VecDestroy(&v3)); 426 PetscCall(VecDestroy(&v4)); 427 PetscCall(VecDestroy(&x)); 428 } 429 430 /* Destroy objects */ 431 PetscCall(VecDestroy(&v1)); 432 PetscCall(VecDestroy(&v2)); 433 PetscCall(PetscRandomDestroy(&rdm)); 434 PetscCall(PetscFree(idxn)); 435 436 PetscCall(MatDestroy(&A_save)); 437 PetscCall(MatDestroy(&B)); 438 439 PetscPreLoadEnd(); 440 PetscCall(PetscFinalize()); 441 return 0; 442 } 443 444 /*TEST 445 446 test: 447 suffix: 2_mattransposematmult_matmatmult 448 nsize: 3 449 requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES) 450 args: -f0 ${DATAFILESPATH}/matrices/medium -f1 ${DATAFILESPATH}/matrices/medium -mattransposematmult_via at*b> ex94_2.tmp 2>&1 451 452 test: 453 suffix: 2_mattransposematmult_scalable 454 nsize: 3 455 requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES) 456 args: -f0 ${DATAFILESPATH}/matrices/medium -f1 ${DATAFILESPATH}/matrices/medium -mattransposematmult_via scalable> ex94_2.tmp 2>&1 457 output_file: output/ex94_1.out 458 459 test: 460 suffix: axpy_mpiaij 461 requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES) 462 nsize: 8 463 args: -f0 ${DATAFILESPATH}/matrices/poisson_2d5p -f1 ${DATAFILESPATH}/matrices/poisson_2d13p -test_MatAXPY 464 output_file: output/ex94_1.out 465 466 test: 467 suffix: axpy_mpibaij 468 requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES) 469 nsize: 8 470 args: -f0 ${DATAFILESPATH}/matrices/poisson_2d5p -f1 ${DATAFILESPATH}/matrices/poisson_2d13p -test_MatAXPY -mat_type baij 471 output_file: output/ex94_1.out 472 473 test: 474 suffix: axpy_mpisbaij 475 requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES) 476 nsize: 8 477 args: -f0 ${DATAFILESPATH}/matrices/poisson_2d5p -f1 ${DATAFILESPATH}/matrices/poisson_2d13p -test_MatAXPY -mat_type sbaij 478 output_file: output/ex94_1.out 479 480 test: 481 suffix: matmatmult 482 requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES) 483 args: -f0 ${DATAFILESPATH}/matrices/arco1 -f1 ${DATAFILESPATH}/matrices/arco1 -viewer_binary_skip_info 484 output_file: output/ex94_1.out 485 486 test: 487 suffix: matmatmult_2 488 requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES) 489 args: -f0 ${DATAFILESPATH}/matrices/arco1 -f1 ${DATAFILESPATH}/matrices/arco1 -mat_type mpiaij -viewer_binary_skip_info 490 output_file: output/ex94_1.out 491 492 test: 493 suffix: matmatmult_scalable 494 nsize: 4 495 requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES) 496 args: -f0 ${DATAFILESPATH}/matrices/arco1 -f1 ${DATAFILESPATH}/matrices/arco1 -matmatmult_via scalable 497 output_file: output/ex94_1.out 498 499 test: 500 suffix: ptap 501 nsize: 3 502 requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES) 503 args: -f0 ${DATAFILESPATH}/matrices/medium -f1 ${DATAFILESPATH}/matrices/medium -matptap_via scalable 504 output_file: output/ex94_1.out 505 506 test: 507 suffix: rap 508 nsize: 3 509 requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES) 510 args: -f0 ${DATAFILESPATH}/matrices/medium -f1 ${DATAFILESPATH}/matrices/medium 511 output_file: output/ex94_1.out 512 513 test: 514 suffix: scalable0 515 requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES) 516 args: -f0 ${DATAFILESPATH}/matrices/arco1 -f1 ${DATAFILESPATH}/matrices/arco1 -viewer_binary_skip_info 517 output_file: output/ex94_1.out 518 519 test: 520 suffix: scalable1 521 requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES) 522 args: -f0 ${DATAFILESPATH}/matrices/arco1 -f1 ${DATAFILESPATH}/matrices/arco1 -viewer_binary_skip_info -matptap_via scalable 523 output_file: output/ex94_1.out 524 525 test: 526 suffix: view 527 nsize: 2 528 requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES) 529 args: -f0 ${DATAFILESPATH}/matrices/tiny -f1 ${DATAFILESPATH}/matrices/tiny -viewer_binary_skip_info -matops_view 530 output_file: output/ex94_2.out 531 532 TEST*/ 533