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