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