1 2 static char help[] = "Tests MatMatSolve() and MatMatTransposeSolve() for computing inv(A) with MUMPS.\n\ 3 Example: mpiexec -n <np> ./ex214 -displ \n\n"; 4 5 #include <petscmat.h> 6 7 int main(int argc,char **args) 8 { 9 PetscMPIInt size,rank; 10 #if defined(PETSC_HAVE_MUMPS) 11 Mat A,RHS,C,F,X,AX,spRHST; 12 PetscInt m,n,nrhs,M,N,i,Istart,Iend,Ii,j,J,test; 13 PetscScalar v; 14 PetscReal norm,tol=PETSC_SQRT_MACHINE_EPSILON; 15 PetscRandom rand; 16 PetscBool displ=PETSC_FALSE; 17 char solver[256]; 18 #endif 19 20 PetscFunctionBeginUser; 21 PetscCall(PetscInitialize(&argc,&args,(char*)0,help)); 22 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 23 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank)); 24 25 #if !defined(PETSC_HAVE_MUMPS) 26 if (rank == 0) PetscCall(PetscPrintf(PETSC_COMM_SELF,"This example requires MUMPS, exit...\n")); 27 PetscCall(PetscFinalize()); 28 return 0; 29 #else 30 31 PetscCall(PetscOptionsGetBool(NULL,NULL,"-displ",&displ,NULL)); 32 33 /* Create matrix A */ 34 m = 4; n = 4; 35 PetscCall(PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL)); 36 PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL)); 37 38 PetscCall(MatCreate(PETSC_COMM_WORLD,&A)); 39 PetscCall(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m*n,m*n)); 40 PetscCall(MatSetFromOptions(A)); 41 PetscCall(MatMPIAIJSetPreallocation(A,5,NULL,5,NULL)); 42 PetscCall(MatSeqAIJSetPreallocation(A,5,NULL)); 43 44 PetscCall(MatGetOwnershipRange(A,&Istart,&Iend)); 45 for (Ii=Istart; Ii<Iend; Ii++) { 46 v = -1.0; i = Ii/n; j = Ii - i*n; 47 if (i>0) {J = Ii - n; PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES));} 48 if (i<m-1) {J = Ii + n; PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES));} 49 if (j>0) {J = Ii - 1; PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES));} 50 if (j<n-1) {J = Ii + 1; PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES));} 51 v = 4.0; PetscCall(MatSetValues(A,1,&Ii,1,&Ii,&v,ADD_VALUES)); 52 } 53 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 54 PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 55 56 PetscCall(MatGetLocalSize(A,&m,&n)); 57 PetscCall(MatGetSize(A,&M,&N)); 58 PetscCheck(m == n,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ, "This example is not intended for rectangular matrices (%" PetscInt_FMT ", %" PetscInt_FMT ")", m, n); 59 60 /* Create dense matrix C and X; C holds true solution with identical columns */ 61 nrhs = N; 62 PetscCall(PetscOptionsGetInt(NULL,NULL,"-nrhs",&nrhs,NULL)); 63 PetscCall(MatCreate(PETSC_COMM_WORLD,&C)); 64 PetscCall(MatSetSizes(C,m,PETSC_DECIDE,PETSC_DECIDE,nrhs)); 65 PetscCall(MatSetType(C,MATDENSE)); 66 PetscCall(MatSetFromOptions(C)); 67 PetscCall(MatSetUp(C)); 68 69 PetscCall(PetscRandomCreate(PETSC_COMM_WORLD,&rand)); 70 PetscCall(PetscRandomSetFromOptions(rand)); 71 PetscCall(MatSetRandom(C,rand)); 72 PetscCall(MatDuplicate(C,MAT_DO_NOT_COPY_VALUES,&X)); 73 74 PetscCall(PetscStrcpy(solver,MATSOLVERMUMPS)); 75 if (rank == 0 && displ) PetscCall(PetscPrintf(PETSC_COMM_SELF,"Solving with %s: nrhs %" PetscInt_FMT ", size mat %" PetscInt_FMT " x %" PetscInt_FMT "\n",solver,nrhs,M,N)); 76 77 for (test=0; test<2; test++) { 78 if (test == 0) { 79 /* Test LU Factorization */ 80 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"test LU factorization\n")); 81 PetscCall(MatGetFactor(A,solver,MAT_FACTOR_LU,&F)); 82 PetscCall(MatLUFactorSymbolic(F,A,NULL,NULL,NULL)); 83 PetscCall(MatLUFactorNumeric(F,A,NULL)); 84 } else { 85 /* Test Cholesky Factorization */ 86 PetscBool flg; 87 PetscCall(MatIsSymmetric(A,0.0,&flg)); 88 PetscCheck(flg,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"A must be symmetric for Cholesky factorization"); 89 90 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"test Cholesky factorization\n")); 91 PetscCall(MatGetFactor(A,solver,MAT_FACTOR_CHOLESKY,&F)); 92 PetscCall(MatCholeskyFactorSymbolic(F,A,NULL,NULL)); 93 PetscCall(MatCholeskyFactorNumeric(F,A,NULL)); 94 } 95 96 /* (1) Test MatMatSolve(): dense RHS = A*C, C: true solutions */ 97 /* ---------------------------------------------------------- */ 98 PetscCall(MatMatMult(A,C,MAT_INITIAL_MATRIX,2.0,&RHS)); 99 PetscCall(MatMatSolve(F,RHS,X)); 100 /* Check the error */ 101 PetscCall(MatAXPY(X,-1.0,C,SAME_NONZERO_PATTERN)); 102 PetscCall(MatNorm(X,NORM_FROBENIUS,&norm)); 103 if (norm > tol) { 104 PetscCall(PetscPrintf(PETSC_COMM_SELF,"(1) MatMatSolve: Norm of error %g\n",norm)); 105 } 106 107 /* Test X=RHS */ 108 PetscCall(MatMatSolve(F,RHS,RHS)); 109 /* Check the error */ 110 PetscCall(MatAXPY(RHS,-1.0,C,SAME_NONZERO_PATTERN)); 111 PetscCall(MatNorm(RHS,NORM_FROBENIUS,&norm)); 112 if (norm > tol) { 113 PetscCall(PetscPrintf(PETSC_COMM_SELF,"(1.1) MatMatSolve: Norm of error %g\n",norm)); 114 } 115 116 /* (2) Test MatMatSolve() for inv(A) with dense RHS: 117 RHS = [e[0],...,e[nrhs-1]], dense X holds first nrhs columns of inv(A) */ 118 /* -------------------------------------------------------------------- */ 119 PetscCall(MatZeroEntries(RHS)); 120 for (i=0; i<nrhs; i++) { 121 v = 1.0; 122 PetscCall(MatSetValues(RHS,1,&i,1,&i,&v,INSERT_VALUES)); 123 } 124 PetscCall(MatAssemblyBegin(RHS,MAT_FINAL_ASSEMBLY)); 125 PetscCall(MatAssemblyEnd(RHS,MAT_FINAL_ASSEMBLY)); 126 127 PetscCall(MatMatSolve(F,RHS,X)); 128 if (displ) { 129 if (rank == 0) PetscCall(PetscPrintf(PETSC_COMM_SELF," \n(2) first %" PetscInt_FMT " columns of inv(A) with dense RHS:\n",nrhs)); 130 PetscCall(MatView(X,PETSC_VIEWER_STDOUT_WORLD)); 131 } 132 133 /* Check the residual */ 134 PetscCall(MatMatMult(A,X,MAT_INITIAL_MATRIX,2.0,&AX)); 135 PetscCall(MatAXPY(AX,-1.0,RHS,SAME_NONZERO_PATTERN)); 136 PetscCall(MatNorm(AX,NORM_INFINITY,&norm)); 137 if (norm > tol) { 138 PetscCall(PetscPrintf(PETSC_COMM_SELF,"(2) MatMatSolve: Norm of residual %g\n",norm)); 139 } 140 PetscCall(MatZeroEntries(X)); 141 142 /* (3) Test MatMatTransposeSolve() for inv(A) with sparse RHS stored in the host: 143 spRHST = [e[0],...,e[nrhs-1]]^T, dense X holds first nrhs columns of inv(A) */ 144 /* --------------------------------------------------------------------------- */ 145 /* Create spRHST: PETSc does not support compressed column format which is required by MUMPS for sparse RHS matrix, 146 thus user must create spRHST=spRHS^T and call MatMatTransposeSolve() */ 147 PetscCall(MatCreate(PETSC_COMM_WORLD,&spRHST)); 148 if (rank == 0) { 149 /* MUMPS requirs RHS be centralized on the host! */ 150 PetscCall(MatSetSizes(spRHST,nrhs,M,PETSC_DECIDE,PETSC_DECIDE)); 151 } else { 152 PetscCall(MatSetSizes(spRHST,0,0,PETSC_DECIDE,PETSC_DECIDE)); 153 } 154 PetscCall(MatSetType(spRHST,MATAIJ)); 155 PetscCall(MatSetFromOptions(spRHST)); 156 PetscCall(MatSetUp(spRHST)); 157 if (rank == 0) { 158 v = 1.0; 159 for (i=0; i<nrhs; i++) { 160 PetscCall(MatSetValues(spRHST,1,&i,1,&i,&v,INSERT_VALUES)); 161 } 162 } 163 PetscCall(MatAssemblyBegin(spRHST,MAT_FINAL_ASSEMBLY)); 164 PetscCall(MatAssemblyEnd(spRHST,MAT_FINAL_ASSEMBLY)); 165 166 PetscCall(MatMatTransposeSolve(F,spRHST,X)); 167 168 if (displ) { 169 if (rank == 0) PetscCall(PetscPrintf(PETSC_COMM_SELF," \n(3) first %" PetscInt_FMT " columns of inv(A) with sparse RHS:\n",nrhs)); 170 PetscCall(MatView(X,PETSC_VIEWER_STDOUT_WORLD)); 171 } 172 173 /* Check the residual: R = A*X - RHS */ 174 PetscCall(MatMatMult(A,X,MAT_REUSE_MATRIX,2.0,&AX)); 175 176 PetscCall(MatAXPY(AX,-1.0,RHS,SAME_NONZERO_PATTERN)); 177 PetscCall(MatNorm(AX,NORM_INFINITY,&norm)); 178 if (norm > tol) { 179 PetscCall(PetscPrintf(PETSC_COMM_SELF,"(3) MatMatSolve: Norm of residual %g\n",norm)); 180 } 181 182 /* (4) Test MatMatSolve() for inv(A) with selected entries: 183 input: spRHS gives selected indices; output: spRHS holds selected entries of inv(A) */ 184 /* --------------------------------------------------------------------------------- */ 185 if (nrhs == N) { /* mumps requires nrhs = n */ 186 /* Create spRHS on proc[0] */ 187 Mat spRHS = NULL; 188 189 /* Create spRHS = spRHST^T. Two matrices share internal matrix data structure */ 190 PetscCall(MatCreateTranspose(spRHST,&spRHS)); 191 PetscCall(MatMumpsGetInverse(F,spRHS)); 192 PetscCall(MatDestroy(&spRHS)); 193 194 PetscCall(MatMumpsGetInverseTranspose(F,spRHST)); 195 if (displ) { 196 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nSelected entries of inv(A^T):\n")); 197 PetscCall(MatView(spRHST,PETSC_VIEWER_STDOUT_WORLD)); 198 } 199 PetscCall(MatDestroy(&spRHS)); 200 } 201 PetscCall(MatDestroy(&AX)); 202 PetscCall(MatDestroy(&F)); 203 PetscCall(MatDestroy(&RHS)); 204 PetscCall(MatDestroy(&spRHST)); 205 } 206 207 /* Free data structures */ 208 PetscCall(MatDestroy(&A)); 209 PetscCall(MatDestroy(&C)); 210 PetscCall(MatDestroy(&X)); 211 PetscCall(PetscRandomDestroy(&rand)); 212 PetscCall(PetscFinalize()); 213 return 0; 214 #endif 215 } 216 217 /*TEST 218 219 test: 220 requires: mumps double !complex 221 222 test: 223 suffix: 2 224 requires: mumps double !complex 225 nsize: 2 226 227 TEST*/ 228