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