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