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