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