static char help[] = "Test Matrix products for AIJ matrices\n\ Input arguments are:\n\ -fA -fB -fC : file to load\n\n"; /* Example of usage: ./ex62 -fA -fB mpiexec -n 3 ./ex62 -fA medium -fB medium */ #include /* B = A - B norm = norm(B) */ PetscErrorCode MatNormDifference(Mat A,Mat B,PetscReal *norm) { PetscFunctionBegin; PetscCall(MatAXPY(B,-1.0,A,DIFFERENT_NONZERO_PATTERN)); PetscCall(MatNorm(B,NORM_FROBENIUS,norm)); PetscFunctionReturn(0); } int main(int argc,char **args) { Mat A,A_save,B,C,P,C1,R; PetscViewer viewer; PetscMPIInt size,rank; PetscInt i,j,*idxn,PM,PN = PETSC_DECIDE,rstart,rend; PetscReal norm; PetscRandom rdm; char file[2][PETSC_MAX_PATH_LEN] = { "", ""}; PetscScalar *a,rval,alpha; PetscBool Test_MatMatMult=PETSC_TRUE,Test_MatTrMat=PETSC_TRUE,Test_MatMatTr=PETSC_TRUE; PetscBool Test_MatPtAP=PETSC_TRUE,Test_MatRARt=PETSC_TRUE,flg,seqaij,flgA,flgB; MatInfo info; PetscInt nzp = 5; /* num of nonzeros in each row of P */ MatType mattype; const char *deft = MATAIJ; char A_mattype[256], B_mattype[256]; PetscInt mcheck = 10; PetscCall(PetscInitialize(&argc,&args,(char*)0,help)); PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank)); /* Load the matrices A_save and B */ PetscOptionsBegin(PETSC_COMM_WORLD,"","",""); PetscCall(PetscOptionsBool("-test_rart","Test MatRARt","",Test_MatRARt,&Test_MatRARt,NULL)); PetscCall(PetscOptionsInt("-PN","Number of columns of P","",PN,&PN,NULL)); PetscCall(PetscOptionsInt("-mcheck","Number of matmult checks","",mcheck,&mcheck,NULL)); PetscCall(PetscOptionsString("-fA","Path for matrix A","",file[0],file[0],sizeof(file[0]),&flg)); PetscCheck(flg,PETSC_COMM_WORLD,PETSC_ERR_USER_INPUT,"Must indicate a file name for matrix A with the -fA option."); PetscCall(PetscOptionsString("-fB","Path for matrix B","",file[1],file[1],sizeof(file[1]),&flg)); PetscCall(PetscOptionsFList("-A_mat_type","Matrix type","MatSetType",MatList,deft,A_mattype,256,&flgA)); PetscCall(PetscOptionsFList("-B_mat_type","Matrix type","MatSetType",MatList,deft,B_mattype,256,&flgB)); PetscOptionsEnd(); PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[0],FILE_MODE_READ,&viewer)); PetscCall(MatCreate(PETSC_COMM_WORLD,&A_save)); PetscCall(MatLoad(A_save,viewer)); PetscCall(PetscViewerDestroy(&viewer)); if (flg) { PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[1],FILE_MODE_READ,&viewer)); PetscCall(MatCreate(PETSC_COMM_WORLD,&B)); PetscCall(MatLoad(B,viewer)); PetscCall(PetscViewerDestroy(&viewer)); } else { PetscCall(PetscObjectReference((PetscObject)A_save)); B = A_save; } if (flgA) { PetscCall(MatConvert(A_save,A_mattype,MAT_INPLACE_MATRIX,&A_save)); } if (flgB) { PetscCall(MatConvert(B,B_mattype,MAT_INPLACE_MATRIX,&B)); } PetscCall(MatSetFromOptions(A_save)); PetscCall(MatSetFromOptions(B)); PetscCall(MatGetType(B,&mattype)); PetscCall(PetscMalloc(nzp*(sizeof(PetscInt)+sizeof(PetscScalar)),&idxn)); a = (PetscScalar*)(idxn + nzp); PetscCall(PetscRandomCreate(PETSC_COMM_WORLD,&rdm)); PetscCall(PetscRandomSetFromOptions(rdm)); /* 1) MatMatMult() */ /* ----------------*/ if (Test_MatMatMult) { PetscCall(MatDuplicate(A_save,MAT_COPY_VALUES,&A)); /* (1.1) Test developer API */ PetscCall(MatProductCreate(A,B,NULL,&C)); PetscCall(MatSetOptionsPrefix(C,"AB_")); PetscCall(MatProductSetType(C,MATPRODUCT_AB)); PetscCall(MatProductSetAlgorithm(C,MATPRODUCTALGORITHMDEFAULT)); PetscCall(MatProductSetFill(C,PETSC_DEFAULT)); PetscCall(MatProductSetFromOptions(C)); /* we can inquire about MATOP_PRODUCTSYMBOLIC even if the destination matrix type has not been set yet */ PetscCall(MatHasOperation(C,MATOP_PRODUCTSYMBOLIC,&flg)); PetscCall(MatProductSymbolic(C)); PetscCall(MatProductNumeric(C)); PetscCall(MatMatMultEqual(A,B,C,mcheck,&flg)); PetscCheck(flg,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error in C=A*B"); /* Test reuse symbolic C */ alpha = 0.9; PetscCall(MatScale(A,alpha)); PetscCall(MatProductNumeric(C)); PetscCall(MatMatMultEqual(A,B,C,mcheck,&flg)); PetscCheck(flg,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error in C=A*B"); PetscCall(MatDestroy(&C)); /* (1.2) Test user driver */ PetscCall(MatMatMult(A,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&C)); /* Test MAT_REUSE_MATRIX - reuse symbolic C */ alpha = 1.0; for (i=0; i<2; i++) { alpha -= 0.1; PetscCall(MatScale(A,alpha)); PetscCall(MatMatMult(A,B,MAT_REUSE_MATRIX,PETSC_DEFAULT,&C)); } PetscCall(MatMatMultEqual(A,B,C,mcheck,&flg)); PetscCheck(flg,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error: MatMatMult()"); PetscCall(MatDestroy(&A)); /* Test MatProductClear() */ PetscCall(MatProductClear(C)); PetscCall(MatDestroy(&C)); /* Test MatMatMult() for dense and aij matrices */ PetscCall(PetscObjectTypeCompareAny((PetscObject)A,&flg,MATSEQAIJ,MATMPIAIJ,"")); if (flg) { PetscCall(MatConvert(A_save,MATDENSE,MAT_INITIAL_MATRIX,&A)); PetscCall(MatMatMult(A,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&C)); PetscCall(MatDestroy(&C)); PetscCall(MatDestroy(&A)); } } /* Create P and R = P^T */ /* --------------------- */ PetscCall(MatGetSize(B,&PM,NULL)); if (PN < 0) PN = PM/2; PetscCall(MatCreate(PETSC_COMM_WORLD,&P)); PetscCall(MatSetSizes(P,PETSC_DECIDE,PETSC_DECIDE,PM,PN)); PetscCall(MatSetType(P,MATAIJ)); PetscCall(MatSeqAIJSetPreallocation(P,nzp,NULL)); PetscCall(MatMPIAIJSetPreallocation(P,nzp,NULL,nzp,NULL)); PetscCall(MatGetOwnershipRange(P,&rstart,&rend)); for (i=0; i