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