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