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 PetscErrorCode ierr; 11 PetscInt m=5,n,p,i,j,nrows,ncols; 12 PetscScalar *v,*barray,rval; 13 PetscReal norm,tol=1.e5*PETSC_MACHINE_EPSILON; 14 PetscMPIInt size,rank; 15 PetscRandom rand; 16 const PetscInt *rows,*cols; 17 IS isrows,iscols; 18 PetscBool mats_view=PETSC_FALSE; 19 20 ierr = PetscInitialize(&argc,&argv,(char*) 0,help);if (ierr) return ierr; 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 ScaLAPACK matrix A\n")); 37 CHKERRQ(MatCreate(PETSC_COMM_WORLD,&A)); 38 CHKERRQ(MatSetSizes(A,m,n,PETSC_DECIDE,PETSC_DECIDE)); 39 CHKERRQ(MatSetType(A,MATSCALAPACK)); 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,MATSCALAPACK)); 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(MatDuplicate(B,MAT_DO_NOT_COPY_VALUES,&X)); 122 123 /* Cholesky factorization */ 124 /*------------------------*/ 125 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD," Create ScaLAPACK matrix Aher\n")); 126 CHKERRQ(MatHermitianTranspose(A,MAT_INITIAL_MATRIX,&Aher)); 127 CHKERRQ(MatAXPY(Aher,1.0,A,SAME_NONZERO_PATTERN)); /* Aher = A + A^T */ 128 CHKERRQ(MatShift(Aher,100.0)); /* add 100.0 to diagonals of Aher to make it spd */ 129 if (mats_view) { 130 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD, "Aher:\n")); 131 CHKERRQ(MatView(Aher,PETSC_VIEWER_STDOUT_WORLD)); 132 } 133 134 /* Cholesky factorization */ 135 /*------------------------*/ 136 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD," Test Cholesky Solver \n")); 137 /* In-place Cholesky */ 138 /* Create matrix factor G, with a copy of Aher */ 139 CHKERRQ(MatDuplicate(Aher,MAT_COPY_VALUES,&G)); 140 141 /* G = L * L^T */ 142 CHKERRQ(MatCholeskyFactor(G,0,0)); 143 if (mats_view) { 144 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD, "Cholesky Factor G:\n")); 145 CHKERRQ(MatView(G,PETSC_VIEWER_STDOUT_WORLD)); 146 } 147 148 /* Solve L * L^T x = b and L * L^T * X = B */ 149 CHKERRQ(MatSolve(G,b,x)); 150 CHKERRQ(MatMatSolve(G,B,X)); 151 CHKERRQ(MatDestroy(&G)); 152 153 /* Out-place Cholesky */ 154 CHKERRQ(MatGetFactor(Aher,MATSOLVERSCALAPACK,MAT_FACTOR_CHOLESKY,&G)); 155 CHKERRQ(MatCholeskyFactorSymbolic(G,Aher,0,NULL)); 156 CHKERRQ(MatCholeskyFactorNumeric(G,Aher,NULL)); 157 if (mats_view) { 158 CHKERRQ(MatView(G,PETSC_VIEWER_STDOUT_WORLD)); 159 } 160 CHKERRQ(MatSolve(G,b,x)); 161 CHKERRQ(MatMatSolve(G,B,X)); 162 CHKERRQ(MatDestroy(&G)); 163 164 /* Check norm(Aher*x - b) */ 165 CHKERRQ(VecCreate(PETSC_COMM_WORLD,&c)); 166 CHKERRQ(VecSetSizes(c,m,PETSC_DECIDE)); 167 CHKERRQ(VecSetFromOptions(c)); 168 CHKERRQ(MatMult(Aher,x,c)); 169 CHKERRQ(VecAXPY(c,-1.0,b)); 170 CHKERRQ(VecNorm(c,NORM_1,&norm)); 171 if (norm > tol) { 172 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"Warning: ||Aher*x - b||=%g for Cholesky\n",(double)norm)); 173 } 174 175 /* Check norm(Aher*X - B) */ 176 CHKERRQ(MatMatMult(Aher,X,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&C)); 177 CHKERRQ(MatAXPY(C,-1.0,B,SAME_NONZERO_PATTERN)); 178 CHKERRQ(MatNorm(C,NORM_1,&norm)); 179 if (norm > tol) { 180 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"Warning: ||Aher*X - B||=%g for Cholesky\n",(double)norm)); 181 } 182 183 /* LU factorization */ 184 /*------------------*/ 185 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD," Test LU Solver \n")); 186 /* In-place LU */ 187 /* Create matrix factor F, with a copy of A */ 188 CHKERRQ(MatDuplicate(A,MAT_COPY_VALUES,&F)); 189 /* Create vector d to test MatSolveAdd() */ 190 CHKERRQ(VecDuplicate(x,&d)); 191 CHKERRQ(VecCopy(x,d)); 192 193 /* PF=LU factorization */ 194 CHKERRQ(MatLUFactor(F,0,0,NULL)); 195 196 /* Solve LUX = PB */ 197 CHKERRQ(MatSolveAdd(F,b,d,x)); 198 CHKERRQ(MatMatSolve(F,B,X)); 199 CHKERRQ(MatDestroy(&F)); 200 201 /* Check norm(A*X - B) */ 202 CHKERRQ(VecCreate(PETSC_COMM_WORLD,&e)); 203 CHKERRQ(VecSetSizes(e,m,PETSC_DECIDE)); 204 CHKERRQ(VecSetFromOptions(e)); 205 CHKERRQ(MatMult(A,x,c)); 206 CHKERRQ(MatMult(A,d,e)); 207 CHKERRQ(VecAXPY(c,-1.0,e)); 208 CHKERRQ(VecAXPY(c,-1.0,b)); 209 CHKERRQ(VecNorm(c,NORM_1,&norm)); 210 if (norm > tol) { 211 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"Warning: ||A*x - b||=%g for LU\n",(double)norm)); 212 } 213 /* Reuse product C; replace Aher with A */ 214 CHKERRQ(MatProductReplaceMats(A,NULL,NULL,C)); 215 CHKERRQ(MatMatMult(A,X,MAT_REUSE_MATRIX,PETSC_DEFAULT,&C)); 216 CHKERRQ(MatAXPY(C,-1.0,B,SAME_NONZERO_PATTERN)); 217 CHKERRQ(MatNorm(C,NORM_1,&norm)); 218 if (norm > tol) { 219 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"Warning: ||A*X - B||=%g for LU\n",(double)norm)); 220 } 221 222 /* Out-place LU */ 223 CHKERRQ(MatGetFactor(A,MATSOLVERSCALAPACK,MAT_FACTOR_LU,&F)); 224 CHKERRQ(MatLUFactorSymbolic(F,A,0,0,NULL)); 225 CHKERRQ(MatLUFactorNumeric(F,A,NULL)); 226 if (mats_view) { 227 CHKERRQ(MatView(F,PETSC_VIEWER_STDOUT_WORLD)); 228 } 229 CHKERRQ(MatSolve(F,b,x)); 230 CHKERRQ(MatMatSolve(F,B,X)); 231 CHKERRQ(MatDestroy(&F)); 232 233 /* Free space */ 234 CHKERRQ(MatDestroy(&A)); 235 CHKERRQ(MatDestroy(&Aher)); 236 CHKERRQ(MatDestroy(&B)); 237 CHKERRQ(MatDestroy(&C)); 238 CHKERRQ(MatDestroy(&X)); 239 CHKERRQ(VecDestroy(&b)); 240 CHKERRQ(VecDestroy(&c)); 241 CHKERRQ(VecDestroy(&d)); 242 CHKERRQ(VecDestroy(&e)); 243 CHKERRQ(VecDestroy(&x)); 244 CHKERRQ(PetscRandomDestroy(&rand)); 245 ierr = PetscFinalize(); 246 return ierr; 247 } 248 249 /*TEST 250 251 build: 252 requires: scalapack 253 254 test: 255 nsize: 2 256 output_file: output/ex245.out 257 258 test: 259 suffix: 2 260 nsize: 6 261 output_file: output/ex245.out 262 263 TEST*/ 264