1 2 static char help[] = "Tests MatCreateLRC()\n\n"; 3 4 /*T 5 Concepts: Low rank correction 6 7 Processors: n 8 T*/ 9 10 #include <petscmat.h> 11 12 int main(int argc,char **args) 13 { 14 Vec x,b,c=NULL; 15 Mat A,U,V,LR,X,LRe; 16 PetscInt M = 5, N = 7; 17 PetscBool flg; 18 19 PetscCall(PetscInitialize(&argc,&args,(char*)0,help)); 20 PetscCall(PetscOptionsGetInt(NULL,NULL,"-m",&M,NULL)); 21 PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&N,NULL)); 22 23 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 24 Create the sparse matrix 25 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 26 PetscCall(MatCreate(PETSC_COMM_WORLD,&A)); 27 PetscCall(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,M,N)); 28 PetscCall(MatSetOptionsPrefix(A,"A_")); 29 PetscCall(MatSetFromOptions(A)); 30 PetscCall(MatSeqAIJSetPreallocation(A,5,NULL)); 31 PetscCall(MatMPIAIJSetPreallocation(A,5,NULL,5,NULL)); 32 PetscCall(MatSetUp(A)); 33 PetscCall(MatSetRandom(A,NULL)); 34 35 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 36 Create the dense matrices 37 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 38 PetscCall(MatCreate(PETSC_COMM_WORLD,&U)); 39 PetscCall(MatSetSizes(U,PETSC_DECIDE,PETSC_DECIDE,M,3)); 40 PetscCall(MatSetType(U,MATDENSE)); 41 PetscCall(MatSetOptionsPrefix(U,"U_")); 42 PetscCall(MatSetFromOptions(U)); 43 PetscCall(MatSetUp(U)); 44 PetscCall(MatSetRandom(U,NULL)); 45 46 PetscCall(MatCreate(PETSC_COMM_WORLD,&V)); 47 PetscCall(MatSetSizes(V,PETSC_DECIDE,PETSC_DECIDE,N,3)); 48 PetscCall(MatSetType(V,MATDENSE)); 49 PetscCall(MatSetOptionsPrefix(V,"V_")); 50 PetscCall(MatSetFromOptions(V)); 51 PetscCall(MatSetUp(V)); 52 PetscCall(MatSetRandom(V,NULL)); 53 54 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 55 Create a vector to hold the diagonal of C 56 A sequential vector can be created as well on each process 57 It is user responsibility to ensure the data in the vector 58 is consistent across processors 59 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 60 PetscCall(PetscOptionsHasName(NULL,NULL,"-use_c",&flg)); 61 if (flg) { 62 PetscCall(VecCreateMPI(PETSC_COMM_WORLD,PETSC_DECIDE,3,&c)); 63 PetscCall(VecSetRandom(c,NULL)); 64 } 65 66 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 67 Create low rank correction matrix 68 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 69 PetscCall(PetscOptionsHasName(NULL,NULL,"-low_rank",&flg)); 70 if (flg) { 71 /* create a low-rank matrix, with no A-matrix */ 72 PetscCall(MatCreateLRC(NULL,U,c,V,&LR)); 73 PetscCall(MatDestroy(&A)); 74 } else { 75 PetscCall(MatCreateLRC(A,U,c,V,&LR)); 76 } 77 78 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 79 Create the low rank correction matrix explicitly to check for 80 correctness 81 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 82 PetscCall(MatHermitianTranspose(V,MAT_INITIAL_MATRIX,&X)); 83 PetscCall(MatDiagonalScale(X,c,NULL)); 84 PetscCall(MatMatMult(U,X,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&LRe)); 85 PetscCall(MatDestroy(&X)); 86 if (A) { 87 PetscCall(MatAYPX(LRe,1.0,A,DIFFERENT_NONZERO_PATTERN)); 88 } 89 90 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 91 Create test vectors 92 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 93 PetscCall(MatCreateVecs(LR,&x,&b)); 94 PetscCall(VecSetRandom(x,NULL)); 95 PetscCall(MatMult(LR,x,b)); 96 PetscCall(MatMultTranspose(LR,b,x)); 97 PetscCall(VecDestroy(&x)); 98 PetscCall(VecDestroy(&b)); 99 100 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 101 Check correctness 102 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 103 PetscCall(MatMultEqual(LR,LRe,10,&flg)); 104 if (!flg) PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Error in MatMult\n")); 105 #if !defined(PETSC_USE_COMPLEX) 106 PetscCall(MatMultHermitianTransposeEqual(LR,LRe,10,&flg)); 107 if (!flg) PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Error in MatMultTranspose\n")); 108 #endif 109 110 PetscCall(MatDestroy(&A)); 111 PetscCall(MatDestroy(&LRe)); 112 PetscCall(MatDestroy(&U)); 113 PetscCall(MatDestroy(&V)); 114 PetscCall(VecDestroy(&c)); 115 PetscCall(MatDestroy(&LR)); 116 117 /* 118 Always call PetscFinalize() before exiting a program. This routine 119 - finalizes the PETSc libraries as well as MPI 120 - provides summary and diagnostic information if certain runtime 121 options are chosen (e.g., -log_view). 122 */ 123 PetscCall(PetscFinalize()); 124 return 0; 125 } 126 127 /*TEST 128 129 testset: 130 output_file: output/ex102_1.out 131 nsize: {{1 2}} 132 args: -low_rank {{0 1}} -use_c {{0 1}} 133 test: 134 suffix: standard 135 test: 136 suffix: cuda 137 requires: cuda 138 args: -A_mat_type aijcusparse -U_mat_type densecuda -V_mat_type densecuda 139 140 TEST*/ 141