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 PetscErrorCode ierr; 10 PetscMPIInt size,rank; 11 #if defined(PETSC_HAVE_MUMPS) 12 Mat A,RHS,C,F,X,AX,spRHST; 13 PetscInt m,n,nrhs,M,N,i,Istart,Iend,Ii,j,J,test; 14 PetscScalar v; 15 PetscReal norm,tol=PETSC_SQRT_MACHINE_EPSILON; 16 PetscRandom rand; 17 PetscBool displ=PETSC_FALSE; 18 char solver[256]; 19 #endif 20 21 ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr; 22 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRMPI(ierr); 23 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRMPI(ierr); 24 25 #if !defined(PETSC_HAVE_MUMPS) 26 if (!rank) {ierr = PetscPrintf(PETSC_COMM_SELF,"This example requires MUMPS, exit...\n");CHKERRQ(ierr);} 27 ierr = PetscFinalize(); 28 return ierr; 29 #else 30 31 ierr = PetscOptionsGetBool(NULL,NULL,"-displ",&displ,NULL);CHKERRQ(ierr); 32 33 /* Create matrix A */ 34 m = 4; n = 4; 35 ierr = PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL);CHKERRQ(ierr); 36 ierr = PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);CHKERRQ(ierr); 37 38 ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr); 39 ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m*n,m*n);CHKERRQ(ierr); 40 ierr = MatSetFromOptions(A);CHKERRQ(ierr); 41 ierr = MatMPIAIJSetPreallocation(A,5,NULL,5,NULL);CHKERRQ(ierr); 42 ierr = MatSeqAIJSetPreallocation(A,5,NULL);CHKERRQ(ierr); 43 44 ierr = MatGetOwnershipRange(A,&Istart,&Iend);CHKERRQ(ierr); 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; ierr = MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);CHKERRQ(ierr);} 48 if (i<m-1) {J = Ii + n; ierr = MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);CHKERRQ(ierr);} 49 if (j>0) {J = Ii - 1; ierr = MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);CHKERRQ(ierr);} 50 if (j<n-1) {J = Ii + 1; ierr = MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);CHKERRQ(ierr);} 51 v = 4.0; ierr = MatSetValues(A,1,&Ii,1,&Ii,&v,ADD_VALUES);CHKERRQ(ierr); 52 } 53 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 54 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 55 56 ierr = MatGetLocalSize(A,&m,&n);CHKERRQ(ierr); 57 ierr = MatGetSize(A,&M,&N);CHKERRQ(ierr); 58 if (m != n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ, "This example is not intended for rectangular matrices (%d, %d)", m, n); 59 60 /* Create dense matrix C and X; C holds true solution with identical colums */ 61 nrhs = N; 62 ierr = PetscOptionsGetInt(NULL,NULL,"-nrhs",&nrhs,NULL);CHKERRQ(ierr); 63 ierr = MatCreate(PETSC_COMM_WORLD,&C);CHKERRQ(ierr); 64 ierr = MatSetSizes(C,m,PETSC_DECIDE,PETSC_DECIDE,nrhs);CHKERRQ(ierr); 65 ierr = MatSetType(C,MATDENSE);CHKERRQ(ierr); 66 ierr = MatSetFromOptions(C);CHKERRQ(ierr); 67 ierr = MatSetUp(C);CHKERRQ(ierr); 68 69 ierr = PetscRandomCreate(PETSC_COMM_WORLD,&rand);CHKERRQ(ierr); 70 ierr = PetscRandomSetFromOptions(rand);CHKERRQ(ierr); 71 ierr = MatSetRandom(C,rand);CHKERRQ(ierr); 72 ierr = MatDuplicate(C,MAT_DO_NOT_COPY_VALUES,&X);CHKERRQ(ierr); 73 74 ierr = PetscStrcpy(solver,MATSOLVERMUMPS);CHKERRQ(ierr); 75 if (!rank && displ) {ierr = PetscPrintf(PETSC_COMM_SELF,"Solving with %s: nrhs %D, size mat %D x %D\n",solver,nrhs,M,N);CHKERRQ(ierr);} 76 77 for (test=0; test<2; test++) { 78 if (test == 0) { 79 /* Test LU Factorization */ 80 ierr = PetscPrintf(PETSC_COMM_WORLD,"test LU factorization\n");CHKERRQ(ierr); 81 ierr = MatGetFactor(A,solver,MAT_FACTOR_LU,&F);CHKERRQ(ierr); 82 ierr = MatLUFactorSymbolic(F,A,NULL,NULL,NULL);CHKERRQ(ierr); 83 ierr = MatLUFactorNumeric(F,A,NULL);CHKERRQ(ierr); 84 } else { 85 /* Test Cholesky Factorization */ 86 PetscBool flg; 87 ierr = MatIsSymmetric(A,0.0,&flg);CHKERRQ(ierr); 88 if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"A must be symmetric for Cholesky factorization"); 89 90 ierr = PetscPrintf(PETSC_COMM_WORLD,"test Cholesky factorization\n");CHKERRQ(ierr); 91 ierr = MatGetFactor(A,solver,MAT_FACTOR_CHOLESKY,&F);CHKERRQ(ierr); 92 ierr = MatCholeskyFactorSymbolic(F,A,NULL,NULL);CHKERRQ(ierr); 93 ierr = MatCholeskyFactorNumeric(F,A,NULL);CHKERRQ(ierr); 94 } 95 96 /* (1) Test MatMatSolve(): dense RHS = A*C, C: true solutions */ 97 /* ---------------------------------------------------------- */ 98 ierr = MatMatMult(A,C,MAT_INITIAL_MATRIX,2.0,&RHS);CHKERRQ(ierr); 99 ierr = MatMatSolve(F,RHS,X);CHKERRQ(ierr); 100 /* Check the error */ 101 ierr = MatAXPY(X,-1.0,C,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 102 ierr = MatNorm(X,NORM_FROBENIUS,&norm);CHKERRQ(ierr); 103 if (norm > tol) { 104 ierr = PetscPrintf(PETSC_COMM_SELF,"(1) MatMatSolve: Norm of error %g\n",norm);CHKERRQ(ierr); 105 } 106 107 /* Test X=RHS */ 108 ierr = MatMatSolve(F,RHS,RHS);CHKERRQ(ierr); 109 /* Check the error */ 110 ierr = MatAXPY(RHS,-1.0,C,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 111 ierr = MatNorm(RHS,NORM_FROBENIUS,&norm);CHKERRQ(ierr); 112 if (norm > tol) { 113 ierr = PetscPrintf(PETSC_COMM_SELF,"(1.1) MatMatSolve: Norm of error %g\n",norm);CHKERRQ(ierr); 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 ierr = MatZeroEntries(RHS);CHKERRQ(ierr); 120 for (i=0; i<nrhs; i++) { 121 v = 1.0; 122 ierr = MatSetValues(RHS,1,&i,1,&i,&v,INSERT_VALUES);CHKERRQ(ierr); 123 } 124 ierr = MatAssemblyBegin(RHS,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 125 ierr = MatAssemblyEnd(RHS,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 126 127 ierr = MatMatSolve(F,RHS,X);CHKERRQ(ierr); 128 if (displ) { 129 if (!rank) {ierr = PetscPrintf(PETSC_COMM_SELF," \n(2) first %D columns of inv(A) with dense RHS:\n",nrhs);CHKERRQ(ierr);} 130 ierr = MatView(X,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 131 } 132 133 /* Check the residual */ 134 ierr = MatMatMult(A,X,MAT_INITIAL_MATRIX,2.0,&AX);CHKERRQ(ierr); 135 ierr = MatAXPY(AX,-1.0,RHS,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 136 ierr = MatNorm(AX,NORM_INFINITY,&norm);CHKERRQ(ierr); 137 if (norm > tol) { 138 ierr = PetscPrintf(PETSC_COMM_SELF,"(2) MatMatSolve: Norm of residual %g\n",norm);CHKERRQ(ierr); 139 } 140 ierr = MatZeroEntries(X);CHKERRQ(ierr); 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 ierr = MatCreate(PETSC_COMM_WORLD,&spRHST);CHKERRQ(ierr); 148 if (!rank) { 149 /* MUMPS requirs RHS be centralized on the host! */ 150 ierr = MatSetSizes(spRHST,nrhs,M,PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr); 151 } else { 152 ierr = MatSetSizes(spRHST,0,0,PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr); 153 } 154 ierr = MatSetType(spRHST,MATAIJ);CHKERRQ(ierr); 155 ierr = MatSetFromOptions(spRHST);CHKERRQ(ierr); 156 ierr = MatSetUp(spRHST);CHKERRQ(ierr); 157 if (!rank) { 158 v = 1.0; 159 for (i=0; i<nrhs; i++) { 160 ierr = MatSetValues(spRHST,1,&i,1,&i,&v,INSERT_VALUES);CHKERRQ(ierr); 161 } 162 } 163 ierr = MatAssemblyBegin(spRHST,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 164 ierr = MatAssemblyEnd(spRHST,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 165 166 ierr = MatMatTransposeSolve(F,spRHST,X);CHKERRQ(ierr); 167 168 if (displ) { 169 if (!rank) {ierr = PetscPrintf(PETSC_COMM_SELF," \n(3) first %D columns of inv(A) with sparse RHS:\n",nrhs);CHKERRQ(ierr);} 170 ierr = MatView(X,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 171 } 172 173 /* Check the residual: R = A*X - RHS */ 174 ierr = MatMatMult(A,X,MAT_REUSE_MATRIX,2.0,&AX);CHKERRQ(ierr); 175 176 ierr = MatAXPY(AX,-1.0,RHS,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 177 ierr = MatNorm(AX,NORM_INFINITY,&norm);CHKERRQ(ierr); 178 if (norm > tol) { 179 ierr = PetscPrintf(PETSC_COMM_SELF,"(3) MatMatSolve: Norm of residual %g\n",norm);CHKERRQ(ierr); 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 ierr = MatCreateTranspose(spRHST,&spRHS);CHKERRQ(ierr); 191 ierr = MatMumpsGetInverse(F,spRHS);CHKERRQ(ierr); 192 ierr = MatDestroy(&spRHS);CHKERRQ(ierr); 193 194 ierr = MatMumpsGetInverseTranspose(F,spRHST);CHKERRQ(ierr); 195 if (displ) { 196 ierr = PetscPrintf(PETSC_COMM_WORLD,"\nSelected entries of inv(A^T):\n");CHKERRQ(ierr); 197 ierr = MatView(spRHST,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 198 } 199 ierr = MatDestroy(&spRHS);CHKERRQ(ierr); 200 } 201 ierr = MatDestroy(&AX);CHKERRQ(ierr); 202 ierr = MatDestroy(&F);CHKERRQ(ierr); 203 ierr = MatDestroy(&RHS);CHKERRQ(ierr); 204 ierr = MatDestroy(&spRHST);CHKERRQ(ierr); 205 } 206 207 /* Free data structures */ 208 ierr = MatDestroy(&A);CHKERRQ(ierr); 209 ierr = MatDestroy(&C);CHKERRQ(ierr); 210 ierr = MatDestroy(&X);CHKERRQ(ierr); 211 ierr = PetscRandomDestroy(&rand);CHKERRQ(ierr); 212 ierr = PetscFinalize(); 213 return ierr; 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