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) PetscCall(MatAYPX(LRe,1.0,A,DIFFERENT_NONZERO_PATTERN)); 81 82 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 83 Create test vectors 84 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 85 PetscCall(MatCreateVecs(LR,&x,&b)); 86 PetscCall(VecSetRandom(x,NULL)); 87 PetscCall(MatMult(LR,x,b)); 88 PetscCall(MatMultTranspose(LR,b,x)); 89 PetscCall(VecDestroy(&x)); 90 PetscCall(VecDestroy(&b)); 91 92 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 93 Check correctness 94 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 95 PetscCall(MatMultEqual(LR,LRe,10,&flg)); 96 if (!flg) PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Error in MatMult\n")); 97 #if !defined(PETSC_USE_COMPLEX) 98 PetscCall(MatMultHermitianTransposeEqual(LR,LRe,10,&flg)); 99 if (!flg) PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Error in MatMultTranspose\n")); 100 #endif 101 102 PetscCall(MatDestroy(&A)); 103 PetscCall(MatDestroy(&LRe)); 104 PetscCall(MatDestroy(&U)); 105 PetscCall(MatDestroy(&V)); 106 PetscCall(VecDestroy(&c)); 107 PetscCall(MatDestroy(&LR)); 108 109 /* 110 Always call PetscFinalize() before exiting a program. This routine 111 - finalizes the PETSc libraries as well as MPI 112 - provides summary and diagnostic information if certain runtime 113 options are chosen (e.g., -log_view). 114 */ 115 PetscCall(PetscFinalize()); 116 return 0; 117 } 118 119 /*TEST 120 121 testset: 122 output_file: output/ex102_1.out 123 nsize: {{1 2}} 124 args: -low_rank {{0 1}} -use_c {{0 1}} 125 test: 126 suffix: standard 127 test: 128 suffix: cuda 129 requires: cuda 130 args: -A_mat_type aijcusparse -U_mat_type densecuda -V_mat_type densecuda 131 132 TEST*/ 133