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