1 static char help[] = " * Example code testing SeqDense matrices with an LDA (leading dimension of the user-allocated arrray) larger than M.\n"; 2 3 #include <petscksp.h> 4 5 int main(int argc,char **argv) 6 { 7 KSP solver; 8 PC pc; 9 Mat A,B; 10 Vec X,Y,Z; 11 MatScalar *a; 12 PetscScalar *b,*x,*y,*z; 13 PetscReal nrm; 14 PetscErrorCode ierr,size=8,lda=10, i,j; 15 16 ierr = PetscInitialize(&argc,&argv,0,help);if (ierr) return ierr; 17 /* Create matrix and three vectors: these are all normal */ 18 ierr = PetscMalloc1(lda*size,&b);CHKERRQ(ierr); 19 for (i=0; i<size; i++) { 20 for (j=0; j<size; j++) { 21 b[i+j*lda] = rand(); 22 } 23 } 24 ierr = MatCreate(MPI_COMM_SELF,&A);CHKERRQ(ierr); 25 ierr = MatSetSizes(A,size,size,size,size);CHKERRQ(ierr); 26 ierr = MatSetType(A,MATSEQDENSE);CHKERRQ(ierr); 27 ierr = MatSeqDenseSetPreallocation(A,NULL);CHKERRQ(ierr); 28 29 ierr = MatDenseGetArray(A,&a);CHKERRQ(ierr); 30 for (i=0; i<size; i++) { 31 for (j=0; j<size; j++) { 32 a[i+j*size] = b[i+j*lda]; 33 } 34 } 35 ierr = MatDenseRestoreArray(A,&a);CHKERRQ(ierr); 36 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 37 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 38 39 ierr = MatCreate(MPI_COMM_SELF,&B);CHKERRQ(ierr); 40 ierr = MatSetSizes(B,size,size,size,size);CHKERRQ(ierr); 41 ierr = MatSetType(B,MATSEQDENSE);CHKERRQ(ierr); 42 ierr = MatSeqDenseSetPreallocation(B,b);CHKERRQ(ierr); 43 ierr = MatDenseSetLDA(B,lda);CHKERRQ(ierr); 44 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 45 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 46 47 ierr = PetscMalloc1(size,&x);CHKERRQ(ierr); 48 for (i=0; i<size; i++) x[i] = 1.0; 49 ierr = VecCreateSeqWithArray(MPI_COMM_SELF,1,size,x,&X);CHKERRQ(ierr); 50 ierr = VecAssemblyBegin(X);CHKERRQ(ierr); 51 ierr = VecAssemblyEnd(X);CHKERRQ(ierr); 52 53 ierr = PetscMalloc1(size,&y);CHKERRQ(ierr); 54 ierr = VecCreateSeqWithArray(MPI_COMM_SELF,1,size,y,&Y);CHKERRQ(ierr); 55 ierr = VecAssemblyBegin(Y);CHKERRQ(ierr); 56 ierr = VecAssemblyEnd(Y);CHKERRQ(ierr); 57 58 ierr = PetscMalloc1(size,&z);CHKERRQ(ierr); 59 ierr = VecCreateSeqWithArray(MPI_COMM_SELF,1,size,z,&Z);CHKERRQ(ierr); 60 ierr = VecAssemblyBegin(Z);CHKERRQ(ierr); 61 ierr = VecAssemblyEnd(Z);CHKERRQ(ierr); 62 63 /* 64 * Solve with A and B 65 */ 66 ierr = KSPCreate(MPI_COMM_SELF,&solver);CHKERRQ(ierr); 67 ierr = KSPSetType(solver,KSPPREONLY);CHKERRQ(ierr); 68 ierr = KSPGetPC(solver,&pc);CHKERRQ(ierr); 69 ierr = PCSetType(pc,PCLU);CHKERRQ(ierr); 70 ierr = KSPSetOperators(solver,A,A);CHKERRQ(ierr); 71 ierr = KSPSolve(solver,X,Y);CHKERRQ(ierr); 72 ierr = KSPSetOperators(solver,B,B);CHKERRQ(ierr); 73 ierr = KSPSolve(solver,X,Z);CHKERRQ(ierr); 74 ierr = VecAXPY(Z,-1.0,Y);CHKERRQ(ierr); 75 ierr = VecNorm(Z,NORM_2,&nrm);CHKERRQ(ierr); 76 ierr = PetscPrintf(PETSC_COMM_SELF,"Test1; error norm=%e\n",nrm);CHKERRQ(ierr); 77 78 /* Free spaces */ 79 ierr = PetscFree(b);CHKERRQ(ierr); 80 ierr = PetscFree(x);CHKERRQ(ierr); 81 ierr = PetscFree(y);CHKERRQ(ierr); 82 ierr = PetscFree(z);CHKERRQ(ierr); 83 ierr = VecDestroy(&X);CHKERRQ(ierr); 84 ierr = VecDestroy(&Y);CHKERRQ(ierr); 85 ierr = VecDestroy(&Z);CHKERRQ(ierr); 86 ierr = MatDestroy(&A);CHKERRQ(ierr); 87 ierr = MatDestroy(&B);CHKERRQ(ierr); 88 ierr = KSPDestroy(&solver);CHKERRQ(ierr); 89 90 ierr = PetscFinalize(); 91 return ierr; 92 } 93 94 /*TEST 95 96 test: 97 98 TEST*/ 99