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