1 2 static char help[] = "Tests LU, Cholesky factorization and MatMatSolve() for a ScaLAPACK dense matrix.\n\n"; 3 4 #include <petscmat.h> 5 6 int main(int argc,char **argv) 7 { 8 Mat A,F,B,X,C,Aher,G; 9 Vec b,x,c,d,e; 10 PetscInt m=5,n,p,i,j,nrows,ncols; 11 PetscScalar *v,*barray,rval; 12 PetscReal norm,tol=1.e5*PETSC_MACHINE_EPSILON; 13 PetscMPIInt size,rank; 14 PetscRandom rand; 15 const PetscInt *rows,*cols; 16 IS isrows,iscols; 17 PetscBool mats_view=PETSC_FALSE; 18 19 PetscCall(PetscInitialize(&argc,&argv,(char*) 0,help)); 20 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank)); 21 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 22 23 PetscCall(PetscRandomCreate(PETSC_COMM_WORLD,&rand)); 24 PetscCall(PetscRandomSetFromOptions(rand)); 25 26 /* Get local dimensions of matrices */ 27 PetscCall(PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL)); 28 n = m; 29 PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL)); 30 p = m/2; 31 PetscCall(PetscOptionsGetInt(NULL,NULL,"-p",&p,NULL)); 32 PetscCall(PetscOptionsHasName(NULL,NULL,"-mats_view",&mats_view)); 33 34 /* Create matrix A */ 35 PetscCall(PetscPrintf(PETSC_COMM_WORLD," Create ScaLAPACK matrix A\n")); 36 PetscCall(MatCreate(PETSC_COMM_WORLD,&A)); 37 PetscCall(MatSetSizes(A,m,n,PETSC_DECIDE,PETSC_DECIDE)); 38 PetscCall(MatSetType(A,MATSCALAPACK)); 39 PetscCall(MatSetFromOptions(A)); 40 PetscCall(MatSetUp(A)); 41 /* Set local matrix entries */ 42 PetscCall(MatGetOwnershipIS(A,&isrows,&iscols)); 43 PetscCall(ISGetLocalSize(isrows,&nrows)); 44 PetscCall(ISGetIndices(isrows,&rows)); 45 PetscCall(ISGetLocalSize(iscols,&ncols)); 46 PetscCall(ISGetIndices(iscols,&cols)); 47 PetscCall(PetscMalloc1(nrows*ncols,&v)); 48 for (i=0;i<nrows;i++) { 49 for (j=0;j<ncols;j++) { 50 PetscCall(PetscRandomGetValue(rand,&rval)); 51 v[i*ncols+j] = rval; 52 } 53 } 54 PetscCall(MatSetValues(A,nrows,rows,ncols,cols,v,INSERT_VALUES)); 55 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 56 PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 57 PetscCall(ISRestoreIndices(isrows,&rows)); 58 PetscCall(ISRestoreIndices(iscols,&cols)); 59 PetscCall(ISDestroy(&isrows)); 60 PetscCall(ISDestroy(&iscols)); 61 PetscCall(PetscFree(v)); 62 if (mats_view) { 63 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "A: nrows %" PetscInt_FMT ", m %" PetscInt_FMT "; ncols %" PetscInt_FMT ", n %" PetscInt_FMT "\n",nrows,m,ncols,n)); 64 PetscCall(MatView(A,PETSC_VIEWER_STDOUT_WORLD)); 65 } 66 67 /* Create rhs matrix B */ 68 PetscCall(PetscPrintf(PETSC_COMM_WORLD," Create rhs matrix B\n")); 69 PetscCall(MatCreate(PETSC_COMM_WORLD,&B)); 70 PetscCall(MatSetSizes(B,m,p,PETSC_DECIDE,PETSC_DECIDE)); 71 PetscCall(MatSetType(B,MATSCALAPACK)); 72 PetscCall(MatSetFromOptions(B)); 73 PetscCall(MatSetUp(B)); 74 PetscCall(MatGetOwnershipIS(B,&isrows,&iscols)); 75 PetscCall(ISGetLocalSize(isrows,&nrows)); 76 PetscCall(ISGetIndices(isrows,&rows)); 77 PetscCall(ISGetLocalSize(iscols,&ncols)); 78 PetscCall(ISGetIndices(iscols,&cols)); 79 PetscCall(PetscMalloc1(nrows*ncols,&v)); 80 for (i=0;i<nrows;i++) { 81 for (j=0;j<ncols;j++) { 82 PetscCall(PetscRandomGetValue(rand,&rval)); 83 v[i*ncols+j] = rval; 84 } 85 } 86 PetscCall(MatSetValues(B,nrows,rows,ncols,cols,v,INSERT_VALUES)); 87 PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 88 PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 89 PetscCall(ISRestoreIndices(isrows,&rows)); 90 PetscCall(ISRestoreIndices(iscols,&cols)); 91 PetscCall(ISDestroy(&isrows)); 92 PetscCall(ISDestroy(&iscols)); 93 PetscCall(PetscFree(v)); 94 if (mats_view) { 95 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "B: nrows %" PetscInt_FMT ", m %" PetscInt_FMT "; ncols %" PetscInt_FMT ", p %" PetscInt_FMT "\n",nrows,m,ncols,p)); 96 PetscCall(MatView(B,PETSC_VIEWER_STDOUT_WORLD)); 97 } 98 99 /* Create rhs vector b and solution x (same size as b) */ 100 PetscCall(VecCreate(PETSC_COMM_WORLD,&b)); 101 PetscCall(VecSetSizes(b,m,PETSC_DECIDE)); 102 PetscCall(VecSetFromOptions(b)); 103 PetscCall(VecGetArray(b,&barray)); 104 for (j=0;j<m;j++) { 105 PetscCall(PetscRandomGetValue(rand,&rval)); 106 barray[j] = rval; 107 } 108 PetscCall(VecRestoreArray(b,&barray)); 109 PetscCall(VecAssemblyBegin(b)); 110 PetscCall(VecAssemblyEnd(b)); 111 if (mats_view) { 112 PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d] b: m %" PetscInt_FMT "\n",rank,m)); 113 PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT)); 114 PetscCall(VecView(b,PETSC_VIEWER_STDOUT_WORLD)); 115 } 116 PetscCall(VecDuplicate(b,&x)); 117 118 /* Create matrix X - same size as B */ 119 PetscCall(PetscPrintf(PETSC_COMM_WORLD," Create solution matrix X\n")); 120 PetscCall(MatDuplicate(B,MAT_DO_NOT_COPY_VALUES,&X)); 121 122 /* Cholesky factorization */ 123 /*------------------------*/ 124 PetscCall(PetscPrintf(PETSC_COMM_WORLD," Create ScaLAPACK matrix Aher\n")); 125 PetscCall(MatHermitianTranspose(A,MAT_INITIAL_MATRIX,&Aher)); 126 PetscCall(MatAXPY(Aher,1.0,A,SAME_NONZERO_PATTERN)); /* Aher = A + A^T */ 127 PetscCall(MatShift(Aher,100.0)); /* add 100.0 to diagonals of Aher to make it spd */ 128 if (mats_view) { 129 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Aher:\n")); 130 PetscCall(MatView(Aher,PETSC_VIEWER_STDOUT_WORLD)); 131 } 132 133 /* Cholesky factorization */ 134 /*------------------------*/ 135 PetscCall(PetscPrintf(PETSC_COMM_WORLD," Test Cholesky Solver \n")); 136 /* In-place Cholesky */ 137 /* Create matrix factor G, with a copy of Aher */ 138 PetscCall(MatDuplicate(Aher,MAT_COPY_VALUES,&G)); 139 140 /* G = L * L^T */ 141 PetscCall(MatCholeskyFactor(G,0,0)); 142 if (mats_view) { 143 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Cholesky Factor G:\n")); 144 PetscCall(MatView(G,PETSC_VIEWER_STDOUT_WORLD)); 145 } 146 147 /* Solve L * L^T x = b and L * L^T * X = B */ 148 PetscCall(MatSolve(G,b,x)); 149 PetscCall(MatMatSolve(G,B,X)); 150 PetscCall(MatDestroy(&G)); 151 152 /* Out-place Cholesky */ 153 PetscCall(MatGetFactor(Aher,MATSOLVERSCALAPACK,MAT_FACTOR_CHOLESKY,&G)); 154 PetscCall(MatCholeskyFactorSymbolic(G,Aher,0,NULL)); 155 PetscCall(MatCholeskyFactorNumeric(G,Aher,NULL)); 156 if (mats_view) { 157 PetscCall(MatView(G,PETSC_VIEWER_STDOUT_WORLD)); 158 } 159 PetscCall(MatSolve(G,b,x)); 160 PetscCall(MatMatSolve(G,B,X)); 161 PetscCall(MatDestroy(&G)); 162 163 /* Check norm(Aher*x - b) */ 164 PetscCall(VecCreate(PETSC_COMM_WORLD,&c)); 165 PetscCall(VecSetSizes(c,m,PETSC_DECIDE)); 166 PetscCall(VecSetFromOptions(c)); 167 PetscCall(MatMult(Aher,x,c)); 168 PetscCall(VecAXPY(c,-1.0,b)); 169 PetscCall(VecNorm(c,NORM_1,&norm)); 170 if (norm > tol) { 171 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Warning: ||Aher*x - b||=%g for Cholesky\n",(double)norm)); 172 } 173 174 /* Check norm(Aher*X - B) */ 175 PetscCall(MatMatMult(Aher,X,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&C)); 176 PetscCall(MatAXPY(C,-1.0,B,SAME_NONZERO_PATTERN)); 177 PetscCall(MatNorm(C,NORM_1,&norm)); 178 if (norm > tol) { 179 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Warning: ||Aher*X - B||=%g for Cholesky\n",(double)norm)); 180 } 181 182 /* LU factorization */ 183 /*------------------*/ 184 PetscCall(PetscPrintf(PETSC_COMM_WORLD," Test LU Solver \n")); 185 /* In-place LU */ 186 /* Create matrix factor F, with a copy of A */ 187 PetscCall(MatDuplicate(A,MAT_COPY_VALUES,&F)); 188 /* Create vector d to test MatSolveAdd() */ 189 PetscCall(VecDuplicate(x,&d)); 190 PetscCall(VecCopy(x,d)); 191 192 /* PF=LU factorization */ 193 PetscCall(MatLUFactor(F,0,0,NULL)); 194 195 /* Solve LUX = PB */ 196 PetscCall(MatSolveAdd(F,b,d,x)); 197 PetscCall(MatMatSolve(F,B,X)); 198 PetscCall(MatDestroy(&F)); 199 200 /* Check norm(A*X - B) */ 201 PetscCall(VecCreate(PETSC_COMM_WORLD,&e)); 202 PetscCall(VecSetSizes(e,m,PETSC_DECIDE)); 203 PetscCall(VecSetFromOptions(e)); 204 PetscCall(MatMult(A,x,c)); 205 PetscCall(MatMult(A,d,e)); 206 PetscCall(VecAXPY(c,-1.0,e)); 207 PetscCall(VecAXPY(c,-1.0,b)); 208 PetscCall(VecNorm(c,NORM_1,&norm)); 209 if (norm > tol) { 210 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Warning: ||A*x - b||=%g for LU\n",(double)norm)); 211 } 212 /* Reuse product C; replace Aher with A */ 213 PetscCall(MatProductReplaceMats(A,NULL,NULL,C)); 214 PetscCall(MatMatMult(A,X,MAT_REUSE_MATRIX,PETSC_DEFAULT,&C)); 215 PetscCall(MatAXPY(C,-1.0,B,SAME_NONZERO_PATTERN)); 216 PetscCall(MatNorm(C,NORM_1,&norm)); 217 if (norm > tol) { 218 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Warning: ||A*X - B||=%g for LU\n",(double)norm)); 219 } 220 221 /* Out-place LU */ 222 PetscCall(MatGetFactor(A,MATSOLVERSCALAPACK,MAT_FACTOR_LU,&F)); 223 PetscCall(MatLUFactorSymbolic(F,A,0,0,NULL)); 224 PetscCall(MatLUFactorNumeric(F,A,NULL)); 225 if (mats_view) { 226 PetscCall(MatView(F,PETSC_VIEWER_STDOUT_WORLD)); 227 } 228 PetscCall(MatSolve(F,b,x)); 229 PetscCall(MatMatSolve(F,B,X)); 230 PetscCall(MatDestroy(&F)); 231 232 /* Free space */ 233 PetscCall(MatDestroy(&A)); 234 PetscCall(MatDestroy(&Aher)); 235 PetscCall(MatDestroy(&B)); 236 PetscCall(MatDestroy(&C)); 237 PetscCall(MatDestroy(&X)); 238 PetscCall(VecDestroy(&b)); 239 PetscCall(VecDestroy(&c)); 240 PetscCall(VecDestroy(&d)); 241 PetscCall(VecDestroy(&e)); 242 PetscCall(VecDestroy(&x)); 243 PetscCall(PetscRandomDestroy(&rand)); 244 PetscCall(PetscFinalize()); 245 return 0; 246 } 247 248 /*TEST 249 250 build: 251 requires: scalapack 252 253 test: 254 nsize: 2 255 output_file: output/ex245.out 256 257 test: 258 suffix: 2 259 nsize: 6 260 output_file: output/ex245.out 261 262 TEST*/ 263