1 2 static char help[] = "Tests relaxation for dense matrices.\n\n"; 3 4 #include <petscmat.h> 5 6 int main(int argc,char **args) 7 { 8 Mat C; 9 Vec u,x,b,e; 10 PetscInt i,n = 10,midx[3]; 11 PetscErrorCode ierr; 12 PetscScalar v[3]; 13 PetscReal omega = 1.0,norm; 14 15 ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr; 16 ierr = PetscOptionsGetReal(NULL,NULL,"-omega",&omega,NULL);CHKERRQ(ierr); 17 ierr = PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);CHKERRQ(ierr); 18 19 ierr = MatCreate(PETSC_COMM_SELF,&C);CHKERRQ(ierr); 20 ierr = MatSetSizes(C,n,n,n,n);CHKERRQ(ierr); 21 ierr = MatSetType(C,MATSEQDENSE);CHKERRQ(ierr); 22 ierr = MatSetUp(C);CHKERRQ(ierr); 23 ierr = VecCreateSeq(PETSC_COMM_SELF,n,&b);CHKERRQ(ierr); 24 ierr = VecCreateSeq(PETSC_COMM_SELF,n,&x);CHKERRQ(ierr); 25 ierr = VecCreateSeq(PETSC_COMM_SELF,n,&u);CHKERRQ(ierr); 26 ierr = VecCreateSeq(PETSC_COMM_SELF,n,&e);CHKERRQ(ierr); 27 ierr = VecSet(u,1.0);CHKERRQ(ierr); 28 ierr = VecSet(x,0.0);CHKERRQ(ierr); 29 30 v[0] = -1.; v[1] = 2.; v[2] = -1.; 31 for (i=1; i<n-1; i++) { 32 midx[0] = i-1; midx[1] = i; midx[2] = i+1; 33 ierr = MatSetValues(C,1,&i,3,midx,v,INSERT_VALUES);CHKERRQ(ierr); 34 } 35 i = 0; midx[0] = 0; midx[1] = 1; 36 v[0] = 2.0; v[1] = -1.; 37 ierr = MatSetValues(C,1,&i,2,midx,v,INSERT_VALUES);CHKERRQ(ierr); 38 i = n-1; midx[0] = n-2; midx[1] = n-1; 39 v[0] = -1.0; v[1] = 2.; 40 ierr = MatSetValues(C,1,&i,2,midx,v,INSERT_VALUES);CHKERRQ(ierr); 41 42 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 43 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 44 45 ierr = MatMult(C,u,b);CHKERRQ(ierr); 46 47 for (i=0; i<n; i++) { 48 ierr = MatSOR(C,b,omega,SOR_FORWARD_SWEEP,0.0,1,1,x);CHKERRQ(ierr); 49 ierr = VecWAXPY(e,-1.0,x,u);CHKERRQ(ierr); 50 ierr = VecNorm(e,NORM_2,&norm);CHKERRQ(ierr); 51 ierr = PetscPrintf(PETSC_COMM_SELF,"2-norm of error %g\n",(double)norm);CHKERRQ(ierr); 52 } 53 ierr = MatDestroy(&C);CHKERRQ(ierr); 54 ierr = VecDestroy(&x);CHKERRQ(ierr); 55 ierr = VecDestroy(&b);CHKERRQ(ierr); 56 ierr = VecDestroy(&u);CHKERRQ(ierr); 57 ierr = VecDestroy(&e);CHKERRQ(ierr); 58 ierr = PetscFinalize(); 59 return ierr; 60 } 61 62 /*TEST 63 64 test: 65 66 TEST*/ 67