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