1 2 static char help[] = "Test Matrix products for AIJ matrices\n\ 3 Input arguments are:\n\ 4 -fA <input_file> -fB <input_file> -fC <input_file>: file to load\n\n"; 5 /* Example of usage: 6 ./ex62 -fA <A_binary> -fB <B_binary> 7 mpiexec -n 3 ./ex62 -fA medium -fB medium 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,C,P,C1,R; 29 PetscViewer viewer; 30 PetscErrorCode ierr; 31 PetscMPIInt size,rank; 32 PetscInt i,j,*idxn,M,N,nzp,PN,rstart,rend; 33 PetscReal norm; 34 PetscRandom rdm; 35 char file[2][128]; 36 PetscScalar *a,rval,alpha; 37 PetscBool Test_MatMatMult=PETSC_TRUE,Test_MatTrMat=PETSC_TRUE,Test_MatMatTr=PETSC_TRUE; 38 PetscBool Test_MatPtAP=PETSC_TRUE,Test_MatRARt=PETSC_TRUE,flg,seqaij; 39 MatInfo info; 40 MatType mattype; 41 42 ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr; 43 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr); 44 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr); 45 46 /* Load the matrices A_save and B */ 47 ierr = PetscOptionsGetString(NULL,NULL,"-fA",file[0],sizeof(file[0]),&flg);CHKERRQ(ierr); 48 if (!flg) SETERRQ(PETSC_COMM_WORLD,1,"Must indicate a file name for small matrix A with the -fA option."); 49 ierr = PetscOptionsGetString(NULL,NULL,"-fB",file[1],sizeof(file[1]),&flg);CHKERRQ(ierr); 50 if (!flg) SETERRQ(PETSC_COMM_WORLD,1,"Must indicate a file name for small matrix B with the -fB option."); 51 52 ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[0],FILE_MODE_READ,&viewer);CHKERRQ(ierr); 53 ierr = MatCreate(PETSC_COMM_WORLD,&A_save);CHKERRQ(ierr); 54 ierr = MatSetFromOptions(A_save);CHKERRQ(ierr); 55 ierr = MatLoad(A_save,viewer);CHKERRQ(ierr); 56 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 57 58 ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[1],FILE_MODE_READ,&viewer);CHKERRQ(ierr); 59 ierr = MatCreate(PETSC_COMM_WORLD,&B);CHKERRQ(ierr); 60 ierr = MatSetFromOptions(B);CHKERRQ(ierr); 61 ierr = MatLoad(B,viewer);CHKERRQ(ierr); 62 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 63 64 ierr = MatGetType(B,&mattype);CHKERRQ(ierr); 65 66 ierr = MatGetSize(B,&M,&N);CHKERRQ(ierr); 67 nzp = PetscMax((PetscInt)(0.1*M),5); 68 ierr = PetscMalloc((nzp+1)*(sizeof(PetscInt)+sizeof(PetscScalar)),&idxn);CHKERRQ(ierr); 69 a = (PetscScalar*)(idxn + nzp); 70 71 ierr = PetscRandomCreate(PETSC_COMM_WORLD,&rdm);CHKERRQ(ierr); 72 ierr = PetscRandomSetFromOptions(rdm);CHKERRQ(ierr); 73 74 /* 1) MatMatMult() */ 75 /* ----------------*/ 76 if (Test_MatMatMult) { 77 ierr = MatDuplicate(A_save,MAT_COPY_VALUES,&A);CHKERRQ(ierr); 78 79 /* (1.1) Test developer API */ 80 ierr = MatProductCreate(A,B,NULL,&C);CHKERRQ(ierr); 81 ierr = MatSetOptionsPrefix(C,"AB_");CHKERRQ(ierr); 82 ierr = MatProductSetType(C,MATPRODUCT_AB);CHKERRQ(ierr); 83 ierr = MatProductSetAlgorithm(C,MATPRODUCTALGORITHM_DEFAULT);CHKERRQ(ierr); 84 ierr = MatProductSetFill(C,PETSC_DEFAULT);CHKERRQ(ierr); 85 ierr = MatProductSetFromOptions(C);CHKERRQ(ierr); 86 /* we can inquire about MATOP_PRODUCTSYMBOLIC even if the destination matrix type has not been set yet */ 87 ierr = MatHasOperation(C,MATOP_PRODUCTSYMBOLIC,&flg);CHKERRQ(ierr); 88 ierr = MatProductSymbolic(C);CHKERRQ(ierr); 89 ierr = MatProductNumeric(C);CHKERRQ(ierr); 90 91 /* Test reuse symbolic C */ 92 alpha = 0.9; 93 ierr = MatScale(A,alpha);CHKERRQ(ierr); 94 ierr = MatProductNumeric(C);CHKERRQ(ierr); 95 96 ierr = MatMatMultEqual(A,B,C,10,&flg);CHKERRQ(ierr); 97 if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error in C=A*B"); 98 ierr = MatDestroy(&C);CHKERRQ(ierr); 99 100 /* (1.2) Test user driver */ 101 ierr = MatMatMult(A,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&C);CHKERRQ(ierr); 102 103 /* Test MAT_REUSE_MATRIX - reuse symbolic C */ 104 alpha = 1.0; 105 for (i=0; i<2; i++) { 106 alpha -= 0.1; 107 ierr = MatScale(A,alpha);CHKERRQ(ierr); 108 ierr = MatMatMult(A,B,MAT_REUSE_MATRIX,PETSC_DEFAULT,&C);CHKERRQ(ierr); 109 } 110 ierr = MatMatMultEqual(A,B,C,10,&flg);CHKERRQ(ierr); 111 if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error: MatMatMult()"); 112 ierr = MatDestroy(&A);CHKERRQ(ierr); 113 114 /* Test MatProductClear() */ 115 ierr = MatProductClear(C);CHKERRQ(ierr); 116 ierr = MatDestroy(&C);CHKERRQ(ierr); 117 118 /* Test MatMatMult() for dense and aij matrices */ 119 ierr = MatConvert(A_save,MATDENSE,MAT_INITIAL_MATRIX,&A);CHKERRQ(ierr); 120 ierr = MatMatMult(A,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&C);CHKERRQ(ierr); 121 ierr = MatDestroy(&C);CHKERRQ(ierr); 122 ierr = MatDestroy(&A);CHKERRQ(ierr); 123 124 } 125 126 /* Create P and R = P^T */ 127 /* --------------------- */ 128 PN = M/2; 129 nzp = 5; /* num of nonzeros in each row of P */ 130 ierr = MatCreate(PETSC_COMM_WORLD,&P);CHKERRQ(ierr); 131 ierr = MatSetSizes(P,PETSC_DECIDE,PETSC_DECIDE,M,PN);CHKERRQ(ierr); 132 ierr = MatSetType(P,mattype);CHKERRQ(ierr); 133 ierr = MatSeqAIJSetPreallocation(P,nzp,NULL);CHKERRQ(ierr); 134 ierr = MatMPIAIJSetPreallocation(P,nzp,NULL,nzp,NULL);CHKERRQ(ierr); 135 ierr = MatGetOwnershipRange(P,&rstart,&rend);CHKERRQ(ierr); 136 for (i=0; i<nzp; i++) { 137 ierr = PetscRandomGetValue(rdm,&a[i]);CHKERRQ(ierr); 138 } 139 for (i=rstart; i<rend; i++) { 140 for (j=0; j<nzp; j++) { 141 ierr = PetscRandomGetValue(rdm,&rval);CHKERRQ(ierr); 142 idxn[j] = (PetscInt)(PetscRealPart(rval)*PN); 143 } 144 ierr = MatSetValues(P,1,&i,nzp,idxn,a,ADD_VALUES);CHKERRQ(ierr); 145 } 146 ierr = MatAssemblyBegin(P,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 147 ierr = MatAssemblyEnd(P,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 148 149 ierr = MatTranspose(P,MAT_INITIAL_MATRIX,&R);CHKERRQ(ierr); 150 151 /* 2) MatTransposeMatMult() */ 152 /* ------------------------ */ 153 if (Test_MatTrMat) { 154 /* (2.1) Test developer driver C = P^T*B */ 155 ierr = MatProductCreate(P,B,NULL,&C);CHKERRQ(ierr); 156 ierr = MatSetOptionsPrefix(C,"AtB_");CHKERRQ(ierr); 157 ierr = MatProductSetType(C,MATPRODUCT_AtB);CHKERRQ(ierr); 158 ierr = MatProductSetAlgorithm(C,MATPRODUCTALGORITHM_DEFAULT);CHKERRQ(ierr); 159 ierr = MatProductSetFill(C,PETSC_DEFAULT);CHKERRQ(ierr); 160 ierr = MatProductSetFromOptions(C);CHKERRQ(ierr); 161 ierr = MatProductSymbolic(C);CHKERRQ(ierr); /* equivalent to MatSetUp() */ 162 ierr = MatSetOption(C,MAT_USE_INODES,PETSC_FALSE);CHKERRQ(ierr); /* illustrate how to call MatSetOption() */ 163 ierr = MatProductNumeric(C);CHKERRQ(ierr); 164 ierr = MatProductNumeric(C);CHKERRQ(ierr); /* test reuse symbolic C */ 165 166 ierr = MatTransposeMatMultEqual(P,B,C,10,&flg);CHKERRQ(ierr); 167 if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Error: developer driver C = P^T*B"); 168 ierr = MatDestroy(&C);CHKERRQ(ierr); 169 170 /* (2.2) Test user driver C = P^T*B */ 171 ierr = MatTransposeMatMult(P,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&C);CHKERRQ(ierr); 172 ierr = MatTransposeMatMult(P,B,MAT_REUSE_MATRIX,PETSC_DEFAULT,&C);CHKERRQ(ierr); 173 ierr = MatGetInfo(C,MAT_GLOBAL_SUM,&info);CHKERRQ(ierr); 174 ierr = MatProductClear(C);CHKERRQ(ierr); 175 176 /* Compare P^T*B and R*B */ 177 ierr = MatMatMult(R,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&C1);CHKERRQ(ierr); 178 ierr = MatNormDifference(C,C1,&norm);CHKERRQ(ierr); 179 if (norm > PETSC_SMALL) SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Error in MatTransposeMatMult(): %g",(double)norm); 180 ierr = MatDestroy(&C1);CHKERRQ(ierr); 181 182 /* Test MatDuplicate() of C=P^T*B */ 183 ierr = MatDuplicate(C,MAT_COPY_VALUES,&C1);CHKERRQ(ierr); 184 ierr = MatDestroy(&C1);CHKERRQ(ierr); 185 ierr = MatDestroy(&C);CHKERRQ(ierr); 186 } 187 188 /* 3) MatMatTransposeMult() */ 189 /* ------------------------ */ 190 if (Test_MatMatTr) { 191 /* C = B*R^T */ 192 ierr = PetscObjectTypeCompare((PetscObject)B,MATSEQAIJ,&seqaij);CHKERRQ(ierr); 193 if (size == 1 && seqaij) { 194 ierr = MatMatTransposeMult(B,R,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&C);CHKERRQ(ierr); 195 ierr = MatSetOptionsPrefix(C,"ABt_");CHKERRQ(ierr); /* enable '-ABt_' for matrix C */ 196 ierr = MatGetInfo(C,MAT_GLOBAL_SUM,&info);CHKERRQ(ierr); 197 198 /* Test MAT_REUSE_MATRIX - reuse symbolic C */ 199 ierr = MatMatTransposeMult(B,R,MAT_REUSE_MATRIX,PETSC_DEFAULT,&C);CHKERRQ(ierr); 200 201 /* Check */ 202 ierr = MatMatMult(B,P,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&C1);CHKERRQ(ierr); 203 ierr = MatNormDifference(C,C1,&norm);CHKERRQ(ierr); 204 if (norm > PETSC_SMALL) SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Error in MatMatTransposeMult() %g",(double)norm); 205 ierr = MatDestroy(&C1);CHKERRQ(ierr); 206 ierr = MatDestroy(&C);CHKERRQ(ierr); 207 } 208 } 209 210 /* 4) Test MatPtAP() */ 211 /*-------------------*/ 212 if (Test_MatPtAP) { 213 ierr = MatDuplicate(A_save,MAT_COPY_VALUES,&A);CHKERRQ(ierr); 214 215 /* (4.1) Test developer API */ 216 ierr = MatProductCreate(A,P,NULL,&C);CHKERRQ(ierr); 217 ierr = MatSetOptionsPrefix(C,"PtAP_");CHKERRQ(ierr); 218 ierr = MatProductSetType(C,MATPRODUCT_PtAP);CHKERRQ(ierr); 219 ierr = MatProductSetAlgorithm(C,MATPRODUCTALGORITHM_DEFAULT);CHKERRQ(ierr); 220 ierr = MatProductSetFill(C,PETSC_DEFAULT);CHKERRQ(ierr); 221 ierr = MatProductSetFromOptions(C);CHKERRQ(ierr); 222 ierr = MatProductSymbolic(C);CHKERRQ(ierr); 223 ierr = MatProductNumeric(C);CHKERRQ(ierr); 224 ierr = MatProductNumeric(C);CHKERRQ(ierr); /* reuse symbolic C */ 225 226 ierr = MatPtAPMultEqual(A,P,C,10,&flg);CHKERRQ(ierr); 227 if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Error in MatProduct_PtAP"); 228 ierr = MatDestroy(&C);CHKERRQ(ierr); 229 230 /* (4.2) Test user driver */ 231 ierr = MatPtAP(A,P,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&C);CHKERRQ(ierr); 232 233 /* Test MAT_REUSE_MATRIX - reuse symbolic C */ 234 alpha=1.0; 235 for (i=0; i<2; i++) { 236 alpha -= 0.1; 237 ierr = MatScale(A,alpha);CHKERRQ(ierr); 238 ierr = MatPtAP(A,P,MAT_REUSE_MATRIX,PETSC_DEFAULT,&C);CHKERRQ(ierr); 239 } 240 241 ierr = MatPtAPMultEqual(A,P,C,10,&flg);CHKERRQ(ierr); 242 if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Error in MatPtAP"); 243 244 /* 5) Test MatRARt() */ 245 /* ----------------- */ 246 if (Test_MatRARt) { 247 Mat RARt; 248 ierr = MatTranspose(P,MAT_REUSE_MATRIX,&R);CHKERRQ(ierr); 249 250 /* (5.1) Test developer driver RARt = R*A*Rt */ 251 ierr = MatProductCreate(A,R,NULL,&RARt);CHKERRQ(ierr); 252 ierr = MatSetOptionsPrefix(RARt,"RARt_");CHKERRQ(ierr); 253 ierr = MatProductSetType(RARt,MATPRODUCT_RARt);CHKERRQ(ierr); 254 ierr = MatProductSetAlgorithm(RARt,MATPRODUCTALGORITHM_DEFAULT);CHKERRQ(ierr); 255 ierr = MatProductSetFill(RARt,PETSC_DEFAULT);CHKERRQ(ierr); 256 ierr = MatProductSetFromOptions(RARt);CHKERRQ(ierr); 257 ierr = MatProductSymbolic(RARt);CHKERRQ(ierr); /* equivalent to MatSetUp() */ 258 ierr = MatSetOption(RARt,MAT_USE_INODES,PETSC_FALSE);CHKERRQ(ierr); /* illustrate how to call MatSetOption() */ 259 ierr = MatProductNumeric(RARt);CHKERRQ(ierr); 260 ierr = MatProductNumeric(RARt);CHKERRQ(ierr); /* test reuse symbolic RARt */ 261 ierr = MatDestroy(&RARt);CHKERRQ(ierr); 262 263 /* (2.2) Test user driver RARt = R*A*Rt */ 264 ierr = MatRARt(A,R,MAT_INITIAL_MATRIX,2.0,&RARt);CHKERRQ(ierr); 265 ierr = MatRARt(A,R,MAT_REUSE_MATRIX,2.0,&RARt);CHKERRQ(ierr); 266 267 ierr = MatNormDifference(C,RARt,&norm);CHKERRQ(ierr); 268 if (norm > PETSC_SMALL) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"|PtAP - RARt| = %g",(double)norm); 269 ierr = MatDestroy(&RARt);CHKERRQ(ierr); 270 } 271 272 ierr = MatDestroy(&A);CHKERRQ(ierr); 273 ierr = MatDestroy(&C);CHKERRQ(ierr); 274 } 275 276 /* Destroy objects */ 277 ierr = PetscRandomDestroy(&rdm);CHKERRQ(ierr); 278 ierr = PetscFree(idxn);CHKERRQ(ierr); 279 280 ierr = MatDestroy(&A_save);CHKERRQ(ierr); 281 ierr = MatDestroy(&B);CHKERRQ(ierr); 282 ierr = MatDestroy(&P);CHKERRQ(ierr); 283 ierr = MatDestroy(&R);CHKERRQ(ierr); 284 285 PetscFinalize(); 286 return ierr; 287 } 288 289 /*TEST 290 test: 291 suffix: 1 292 requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES) 293 args: -fA ${DATAFILESPATH}/matrices/medium -fB ${DATAFILESPATH}/matrices/medium 294 output_file: output/ex62_1.out 295 296 test: 297 suffix: 2_ab_scalable 298 requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES) 299 args: -fA ${DATAFILESPATH}/matrices/medium -fB ${DATAFILESPATH}/matrices/medium -AB_matproduct_ab_via scalable -matmatmult_via scalable -AtB_matproduct_atb_via outerproduct -mattransposematmult_via outerproduct 300 output_file: output/ex62_1.out 301 302 test: 303 suffix: 3_ab_scalable_fast 304 requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES) 305 args: -fA ${DATAFILESPATH}/matrices/medium -fB ${DATAFILESPATH}/matrices/medium -AB_matproduct_ab_via scalable_fast -matmatmult_via scalable_fast -matmattransmult_via color 306 output_file: output/ex62_1.out 307 308 test: 309 suffix: 4_ab_heap 310 requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES) 311 args: -fA ${DATAFILESPATH}/matrices/medium -fB ${DATAFILESPATH}/matrices/medium -AB_matproduct_ab_via heap -matmatmult_via heap -PtAP_matproduct_ptap_via rap -matptap_via rap 312 output_file: output/ex62_1.out 313 314 test: 315 suffix: 5_ab_btheap 316 requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES) 317 args: -fA ${DATAFILESPATH}/matrices/medium -fB ${DATAFILESPATH}/matrices/medium -AB_matproduct_ab_via btheap -matmatmult_via btheap -matrart_via r*art 318 output_file: output/ex62_1.out 319 320 test: 321 suffix: 6_ab_llcondensed 322 requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES) 323 args: -fA ${DATAFILESPATH}/matrices/medium -fB ${DATAFILESPATH}/matrices/medium -AB_matproduct_ab_via llcondensed -matmatmult_via llcondensed -matrart_via coloring_rart 324 output_file: output/ex62_1.out 325 326 test: 327 suffix: 7_ab_rowmerge 328 requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES) 329 args: -fA ${DATAFILESPATH}/matrices/medium -fB ${DATAFILESPATH}/matrices/medium -AB_matproduct_ab_via rowmerge -matmatmult_via rowmerge 330 output_file: output/ex62_1.out 331 332 test: 333 suffix: 8_ab_hypre 334 requires: hypre datafilespath !complex double !define(PETSC_USE_64BIT_INDICES) 335 args: -fA ${DATAFILESPATH}/matrices/medium -fB ${DATAFILESPATH}/matrices/medium -AB_matproduct_ab_via hypre -matmatmult_via hypre -PtAP_matproduct_ptap_via hypre -matptap_via hypre 336 output_file: output/ex62_1.out 337 338 test: 339 suffix: 10 340 requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES) 341 nsize: 3 342 args: -fA ${DATAFILESPATH}/matrices/medium -fB ${DATAFILESPATH}/matrices/medium 343 output_file: output/ex62_1.out 344 345 test: 346 suffix: 11_ab_scalable 347 requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES) 348 nsize: 3 349 args: -fA ${DATAFILESPATH}/matrices/medium -fB ${DATAFILESPATH}/matrices/medium -AB_matproduct_ab_via scalable -matmatmult_via scalable -AtB_matproduct_atb_via scalable -mattransposematmult_via scalable 350 output_file: output/ex62_1.out 351 352 test: 353 suffix: 12_ab_seqmpi 354 requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES) 355 nsize: 3 356 args: -fA ${DATAFILESPATH}/matrices/medium -fB ${DATAFILESPATH}/matrices/medium -AB_matproduct_ab_via seqmpi -matmatmult_via seqmpi -AtB_matproduct_atb_via at*b -mattransposematmult_via at*b 357 output_file: output/ex62_1.out 358 359 test: 360 suffix: 13_ab_hypre 361 requires: hypre datafilespath !complex double !define(PETSC_USE_64BIT_INDICES) 362 nsize: 3 363 args: -fA ${DATAFILESPATH}/matrices/medium -fB ${DATAFILESPATH}/matrices/medium -AB_matproduct_ab_via hypre -matmatmult_via hypre -PtAP_matproduct_ptap_via hypre -matptap_via hypre 364 output_file: output/ex62_1.out 365 366 TEST*/ 367