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