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; 16 PetscInt i,j,Ii,J,Istart,Iend,m = 8,n = 7,rstart,rend; 17 PetscErrorCode ierr; 18 PetscBool flg; 19 PetscScalar *u,a; 20 21 ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr; 22 ierr = PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL);CHKERRQ(ierr); 23 ierr = PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);CHKERRQ(ierr); 24 25 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 26 Create the sparse matrix 27 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 28 ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr); 29 ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m*n,m*n);CHKERRQ(ierr); 30 ierr = MatSetFromOptions(A);CHKERRQ(ierr); 31 ierr = MatSetUp(A);CHKERRQ(ierr); 32 ierr = MatGetOwnershipRange(A,&Istart,&Iend);CHKERRQ(ierr); 33 for (Ii=Istart; Ii<Iend; Ii++) { 34 a = -1.0; i = Ii/n; j = Ii - i*n; 35 if (i>0) {J = Ii - n; ierr = MatSetValues(A,1,&Ii,1,&J,&a,INSERT_VALUES);CHKERRQ(ierr);} 36 if (i<m-1) {J = Ii + n; ierr = MatSetValues(A,1,&Ii,1,&J,&a,INSERT_VALUES);CHKERRQ(ierr);} 37 if (j>0) {J = Ii - 1; ierr = MatSetValues(A,1,&Ii,1,&J,&a,INSERT_VALUES);CHKERRQ(ierr);} 38 if (j<n-1) {J = Ii + 1; ierr = MatSetValues(A,1,&Ii,1,&J,&a,INSERT_VALUES);CHKERRQ(ierr);} 39 a = 4.0; ierr = MatSetValues(A,1,&Ii,1,&Ii,&a,INSERT_VALUES);CHKERRQ(ierr); 40 } 41 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 42 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 43 44 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 45 Create the dense matrices 46 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 47 ierr = MatCreate(PETSC_COMM_WORLD,&U);CHKERRQ(ierr); 48 ierr = MatSetSizes(U,PETSC_DECIDE,PETSC_DECIDE,m*n,3);CHKERRQ(ierr); 49 ierr = MatSetType(U,MATDENSE);CHKERRQ(ierr); 50 ierr = MatSetUp(U);CHKERRQ(ierr); 51 ierr = MatGetOwnershipRange(U,&rstart,&rend);CHKERRQ(ierr); 52 ierr = MatDenseGetArray(U,&u);CHKERRQ(ierr); 53 for (i=rstart; i<rend; i++) { 54 u[i-rstart] = (PetscReal)i; 55 u[i+rend-2*rstart] = (PetscReal)1000*i; 56 u[i+2*rend-3*rstart] = (PetscReal)100000*i; 57 } 58 ierr = MatDenseRestoreArray(U,&u);CHKERRQ(ierr); 59 ierr = MatAssemblyBegin(U,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 60 ierr = MatAssemblyEnd(U,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 61 62 ierr = MatCreate(PETSC_COMM_WORLD,&V);CHKERRQ(ierr); 63 ierr = MatSetSizes(V,PETSC_DECIDE,PETSC_DECIDE,m*n,3);CHKERRQ(ierr); 64 ierr = MatSetType(V,MATDENSE);CHKERRQ(ierr); 65 ierr = MatSetUp(V);CHKERRQ(ierr); 66 ierr = MatGetOwnershipRange(U,&rstart,&rend);CHKERRQ(ierr); 67 ierr = MatDenseGetArray(V,&u);CHKERRQ(ierr); 68 for (i=rstart; i<rend; i++) { 69 u[i-rstart] = (PetscReal)i; 70 u[i+rend-2*rstart] = (PetscReal)1.2*i; 71 u[i+2*rend-3*rstart] = (PetscReal)1.67*i+2; 72 } 73 ierr = MatDenseRestoreArray(V,&u);CHKERRQ(ierr); 74 ierr = MatAssemblyBegin(V,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 75 ierr = MatAssemblyEnd(V,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 76 77 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 78 Create a vector to hold the diagonal of C 79 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 80 ierr = PetscOptionsHasName(NULL,NULL,"-use_c",&flg);CHKERRQ(ierr); 81 if (flg) { 82 ierr = VecCreateSeq(PETSC_COMM_SELF,3,&c);CHKERRQ(ierr); 83 ierr = VecGetArray(c,&u);CHKERRQ(ierr); 84 u[0] = 2.0; 85 u[1] = -1.0; 86 u[2] = 1.0; 87 ierr = VecRestoreArray(c,&u);CHKERRQ(ierr); 88 } 89 90 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 91 Create low rank correction matrix 92 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 93 ierr = PetscOptionsHasName(NULL,NULL,"-low_rank",&flg);CHKERRQ(ierr); 94 if (flg) { 95 /* create a low-rank matrix, with no A-matrix */ 96 ierr = MatCreateLRC(NULL,U,c,V,&LR);CHKERRQ(ierr); 97 } else { 98 ierr = MatCreateLRC(A,U,c,V,&LR);CHKERRQ(ierr); 99 } 100 101 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 102 Create test vectors 103 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 104 ierr = VecCreate(PETSC_COMM_WORLD,&x);CHKERRQ(ierr); 105 ierr = VecSetSizes(x,PETSC_DECIDE,m*n);CHKERRQ(ierr); 106 ierr = VecSetFromOptions(x);CHKERRQ(ierr); 107 ierr = VecDuplicate(x,&b);CHKERRQ(ierr); 108 ierr = VecGetOwnershipRange(x,&rstart,&rend);CHKERRQ(ierr); 109 ierr = VecGetArray(x,&u);CHKERRQ(ierr); 110 for (i=rstart; i<rend; i++) u[i-rstart] = (PetscScalar)i; 111 ierr = VecRestoreArray(x,&u);CHKERRQ(ierr); 112 113 ierr = MatMult(LR,x,b);CHKERRQ(ierr); 114 /* 115 View the product if desired 116 */ 117 ierr = PetscOptionsHasName(NULL,NULL,"-view_product",&flg);CHKERRQ(ierr); 118 if (flg) {ierr = VecView(b,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);} 119 120 ierr = VecDestroy(&x);CHKERRQ(ierr); 121 ierr = VecDestroy(&b);CHKERRQ(ierr); 122 /* you can destroy the matrices in any order you like */ 123 ierr = MatDestroy(&A);CHKERRQ(ierr); 124 ierr = MatDestroy(&U);CHKERRQ(ierr); 125 ierr = MatDestroy(&V);CHKERRQ(ierr); 126 ierr = VecDestroy(&c);CHKERRQ(ierr); 127 ierr = MatDestroy(&LR);CHKERRQ(ierr); 128 129 /* 130 Always call PetscFinalize() before exiting a program. This routine 131 - finalizes the PETSc libraries as well as MPI 132 - provides summary and diagnostic information if certain runtime 133 options are chosen (e.g., -log_view). 134 */ 135 ierr = PetscFinalize(); 136 return ierr; 137 } 138 139 140 /*TEST 141 142 test: 143 suffix: 1 144 nsize: 2 145 args: -view_product 146 147 test: 148 suffix: 2 149 nsize: 2 150 args: -low_rank -view_product 151 152 test: 153 suffix: 3 154 nsize: 2 155 args: -use_c -view_product 156 157 TEST*/ 158