1 2 static char help[] = "Tests LU, Cholesky factorization and MatMatSolve() for a sequential dense matrix. \n\ 3 For MATSEQDENSE matrix, the factorization is just a thin wrapper to LAPACK. \n\ 4 For MATSEQDENSECUDA, it uses cusolverDn routines \n\n"; 5 6 #include <petscmat.h> 7 8 int main(int argc,char **argv) 9 { 10 Mat mat,F,RHS,SOLU; 11 MatInfo info; 12 PetscErrorCode ierr; 13 PetscInt n = 10,i,j,rstart,rend,nrhs=2; 14 PetscScalar value = 1.0; 15 Vec x,y,b,ytmp; 16 IS perm; 17 PetscReal norm,tol=PETSC_SMALL; 18 PetscMPIInt size; 19 PetscRandom rand; 20 char solver[64]; 21 PetscBool inplace,full = PETSC_FALSE, ldl = PETSC_TRUE; 22 23 ierr = PetscInitialize(&argc,&argv,(char*) 0,help);if (ierr) return ierr; 24 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr); 25 if (size != 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"This is a uniprocessor example only!"); 26 ierr = PetscStrcpy(solver,"petsc");CHKERRQ(ierr); 27 ierr = PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);CHKERRQ(ierr); 28 ierr = PetscOptionsGetInt(NULL,NULL,"-nrhs",&nrhs,NULL);CHKERRQ(ierr); 29 ierr = PetscOptionsGetBool(NULL,NULL,"-ldl",&ldl,NULL);CHKERRQ(ierr); 30 ierr = PetscOptionsGetBool(NULL,NULL,"-full",&full,NULL);CHKERRQ(ierr); 31 ierr = PetscOptionsGetString(NULL,NULL,"-solver_type",solver,sizeof(solver),NULL);CHKERRQ(ierr); 32 33 /* create multiple vectors RHS and SOLU */ 34 ierr = MatCreate(PETSC_COMM_WORLD,&RHS);CHKERRQ(ierr); 35 ierr = MatSetSizes(RHS,PETSC_DECIDE,PETSC_DECIDE,n,nrhs);CHKERRQ(ierr); 36 ierr = MatSetType(RHS,MATDENSE);CHKERRQ(ierr); 37 ierr = MatSetOptionsPrefix(RHS,"rhs_");CHKERRQ(ierr); 38 ierr = MatSetFromOptions(RHS);CHKERRQ(ierr); 39 ierr = MatSeqDenseSetPreallocation(RHS,NULL);CHKERRQ(ierr); 40 41 ierr = PetscRandomCreate(PETSC_COMM_WORLD,&rand);CHKERRQ(ierr); 42 ierr = PetscRandomSetFromOptions(rand);CHKERRQ(ierr); 43 ierr = MatSetRandom(RHS,rand);CHKERRQ(ierr); 44 45 ierr = MatDuplicate(RHS,MAT_DO_NOT_COPY_VALUES,&SOLU);CHKERRQ(ierr); 46 47 /* create matrix */ 48 ierr = MatCreate(PETSC_COMM_WORLD,&mat);CHKERRQ(ierr); 49 ierr = MatSetSizes(mat,PETSC_DECIDE,PETSC_DECIDE,n,n);CHKERRQ(ierr); 50 ierr = MatSetType(mat,MATDENSE);CHKERRQ(ierr); 51 ierr = MatSetFromOptions(mat);CHKERRQ(ierr); 52 ierr = MatSetUp(mat);CHKERRQ(ierr); 53 ierr = MatGetOwnershipRange(mat,&rstart,&rend);CHKERRQ(ierr); 54 if (!full) { 55 for (i=rstart; i<rend; i++) { 56 value = (PetscReal)i+1; 57 ierr = MatSetValues(mat,1,&i,1,&i,&value,INSERT_VALUES);CHKERRQ(ierr); 58 } 59 ierr = MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 60 ierr = MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 61 } else { 62 Mat T; 63 64 ierr = MatSetRandom(mat,rand);CHKERRQ(ierr); 65 ierr = MatMatTransposeMult(mat,mat,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&T);CHKERRQ(ierr); 66 ierr = MatDestroy(&mat);CHKERRQ(ierr); 67 mat = T; 68 } 69 70 /* create single vectors */ 71 ierr = MatCreateVecs(mat,&x,&b);CHKERRQ(ierr); 72 ierr = VecDuplicate(x,&y);CHKERRQ(ierr); 73 ierr = VecDuplicate(y,&ytmp);CHKERRQ(ierr); 74 ierr = VecSet(x,value);CHKERRQ(ierr); 75 76 /* Only SeqDense* support in-place factorizations and NULL permutations */ 77 ierr = PetscObjectBaseTypeCompare((PetscObject)mat,MATSEQDENSE,&inplace);CHKERRQ(ierr); 78 ierr = MatGetLocalSize(mat,&i,NULL);CHKERRQ(ierr); 79 ierr = MatGetOwnershipRange(mat,&j,NULL);CHKERRQ(ierr); 80 ierr = ISCreateStride(PETSC_COMM_WORLD,i,j,1,&perm);CHKERRQ(ierr); 81 82 ierr = MatGetInfo(mat,MAT_LOCAL,&info);CHKERRQ(ierr); 83 ierr = PetscPrintf(PETSC_COMM_WORLD,"matrix nonzeros = %D, allocated nonzeros = %D\n", 84 (PetscInt)info.nz_used,(PetscInt)info.nz_allocated);CHKERRQ(ierr); 85 ierr = MatMult(mat,x,b);CHKERRQ(ierr); 86 87 /* Cholesky factorization - perm and factinfo are ignored by LAPACK */ 88 /* in-place Cholesky */ 89 if (inplace) { 90 Mat RHS2; 91 92 ierr = MatDuplicate(mat,MAT_COPY_VALUES,&F);CHKERRQ(ierr); 93 if (!ldl) { ierr = MatSetOption(F,MAT_SPD,PETSC_TRUE);CHKERRQ(ierr); } 94 ierr = MatCholeskyFactor(F,perm,0);CHKERRQ(ierr); 95 ierr = MatSolve(F,b,y);CHKERRQ(ierr); 96 ierr = VecAXPY(y,-1.0,x);CHKERRQ(ierr); 97 ierr = VecNorm(y,NORM_2,&norm);CHKERRQ(ierr); 98 if (norm > tol) { 99 ierr = PetscPrintf(PETSC_COMM_WORLD,"Warning: Norm of error for in-place Cholesky %g\n",(double)norm);CHKERRQ(ierr); 100 } 101 102 ierr = MatMatSolve(F,RHS,SOLU);CHKERRQ(ierr); 103 ierr = MatMatMult(mat,SOLU,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&RHS2);CHKERRQ(ierr); 104 ierr = MatAXPY(RHS,-1.0,RHS2,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 105 ierr = MatNorm(RHS,NORM_FROBENIUS,&norm);CHKERRQ(ierr); 106 if (norm > tol) { 107 ierr = PetscPrintf(PETSC_COMM_WORLD,"Error: Norm of residual for in-place Cholesky (MatMatSolve) %g\n",(double)norm);CHKERRQ(ierr); 108 } 109 ierr = MatDestroy(&F);CHKERRQ(ierr); 110 ierr = MatDestroy(&RHS2);CHKERRQ(ierr); 111 } 112 113 /* out-of-place Cholesky */ 114 ierr = MatGetFactor(mat,solver,MAT_FACTOR_CHOLESKY,&F);CHKERRQ(ierr); 115 if (!ldl) { ierr = MatSetOption(F,MAT_SPD,PETSC_TRUE);CHKERRQ(ierr); } 116 ierr = MatCholeskyFactorSymbolic(F,mat,perm,0);CHKERRQ(ierr); 117 ierr = MatCholeskyFactorNumeric(F,mat,0);CHKERRQ(ierr); 118 ierr = MatSolve(F,b,y);CHKERRQ(ierr); 119 ierr = VecAXPY(y,-1.0,x);CHKERRQ(ierr); 120 ierr = VecNorm(y,NORM_2,&norm);CHKERRQ(ierr); 121 if (norm > tol) { 122 ierr = PetscPrintf(PETSC_COMM_WORLD,"Warning: Norm of error for out-of-place Cholesky %g\n",(double)norm);CHKERRQ(ierr); 123 } 124 ierr = MatDestroy(&F);CHKERRQ(ierr); 125 126 /* LU factorization - perms and factinfo are ignored by LAPACK */ 127 i = n-1; 128 ierr = MatZeroRows(mat,1,&i,-1.0,NULL,NULL);CHKERRQ(ierr); 129 ierr = MatMult(mat,x,b);CHKERRQ(ierr); 130 131 /* in-place LU */ 132 if (inplace) { 133 Mat RHS2; 134 135 ierr = MatDuplicate(mat,MAT_COPY_VALUES,&F);CHKERRQ(ierr); 136 ierr = MatLUFactor(F,perm,perm,0);CHKERRQ(ierr); 137 ierr = MatSolve(F,b,y);CHKERRQ(ierr); 138 ierr = VecAXPY(y,-1.0,x);CHKERRQ(ierr); 139 ierr = VecNorm(y,NORM_2,&norm);CHKERRQ(ierr); 140 if (norm > tol) { 141 ierr = PetscPrintf(PETSC_COMM_WORLD,"Warning: Norm of error for in-place LU %g\n",(double)norm);CHKERRQ(ierr); 142 } 143 ierr = MatMatSolve(F,RHS,SOLU);CHKERRQ(ierr); 144 ierr = MatMatMult(mat,SOLU,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&RHS2);CHKERRQ(ierr); 145 ierr = MatAXPY(RHS,-1.0,RHS2,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 146 ierr = MatNorm(RHS,NORM_FROBENIUS,&norm);CHKERRQ(ierr); 147 if (norm > tol) { 148 ierr = PetscPrintf(PETSC_COMM_WORLD,"Error: Norm of residual for in-place LU (MatMatSolve) %g\n",(double)norm);CHKERRQ(ierr); 149 } 150 ierr = MatDestroy(&F);CHKERRQ(ierr); 151 ierr = MatDestroy(&RHS2);CHKERRQ(ierr); 152 } 153 154 /* out-of-place LU */ 155 ierr = MatGetFactor(mat,solver,MAT_FACTOR_LU,&F);CHKERRQ(ierr); 156 ierr = MatLUFactorSymbolic(F,mat,perm,perm,0);CHKERRQ(ierr); 157 ierr = MatLUFactorNumeric(F,mat,0);CHKERRQ(ierr); 158 ierr = MatSolve(F,b,y);CHKERRQ(ierr); 159 ierr = VecAXPY(y,-1.0,x);CHKERRQ(ierr); 160 ierr = VecNorm(y,NORM_2,&norm);CHKERRQ(ierr); 161 if (norm > tol) { 162 ierr = PetscPrintf(PETSC_COMM_WORLD,"Warning: Norm of error for out-of-place LU %g\n",(double)norm);CHKERRQ(ierr); 163 } 164 165 /* free space */ 166 ierr = ISDestroy(&perm);CHKERRQ(ierr); 167 ierr = MatDestroy(&F);CHKERRQ(ierr); 168 ierr = MatDestroy(&mat);CHKERRQ(ierr); 169 ierr = MatDestroy(&RHS);CHKERRQ(ierr); 170 ierr = MatDestroy(&SOLU);CHKERRQ(ierr); 171 ierr = PetscRandomDestroy(&rand);CHKERRQ(ierr); 172 ierr = VecDestroy(&x);CHKERRQ(ierr); 173 ierr = VecDestroy(&b);CHKERRQ(ierr); 174 ierr = VecDestroy(&y);CHKERRQ(ierr); 175 ierr = VecDestroy(&ytmp);CHKERRQ(ierr); 176 ierr = PetscFinalize(); 177 return ierr; 178 } 179 180 181 182 /*TEST 183 184 test: 185 186 test: 187 requires: cuda 188 suffix: seqdensecuda 189 args: -mat_type seqdensecuda -rhs_mat_type seqdensecuda -ldl 0 -solver_type {{petsc cuda}} 190 output_file: output/ex1_1.out 191 192 test: 193 requires: cuda 194 suffix: seqdensecuda_seqaijcusparse 195 args: -mat_type seqaijcusparse -rhs_mat_type seqdensecuda 196 output_file: output/ex1_2.out 197 198 test: 199 requires: cuda viennacl 200 suffix: seqdensecuda_seqaijviennacl 201 args: -mat_type seqaijviennacl -rhs_mat_type seqdensecuda 202 output_file: output/ex1_2.out 203 204 TEST*/ 205