Actual source code: ex3.c
1: /*
2: * Example code testing SeqDense matrices with an LDA (leading dimension
3: * of the user-allocated arrray) larger than M.
4: * This example tests the functionality of MatSolve.
5: */
6: #include <stdlib.h>
7: #include petscmat.h
8: #include petscksp.h
12: int main(int argc,char **argv)
13: {
14: KSP solver; PC pc;
15: Mat A,B;
16: Vec X,Y,Z;
17: PetscScalar *a,*b,*x,*y,*z,one=1,mone=-1;
18: PetscReal nrm;
19: int ierr,size=8,lda=10, i,j;
21: PetscInitialize(&argc,&argv,0,0);
23: /*
24: * Create matrix and three vectors: these are all normal
25: */
26: PetscMalloc(size*size*sizeof(PetscScalar),&a);
27: PetscMalloc(lda*size*sizeof(PetscScalar),&b);
28: for (i=0; i<size; i++) {
29: for (j=0; j<size; j++) {
30: a[i+j*size] = rand(); b[i+j*lda] = a[i+j*size];
31: }
32: }
33: MatCreate(MPI_COMM_SELF,size,size,size,size,&A);
34: MatSetType(A,MATSEQDENSE);
35: MatSeqDenseSetPreallocation(A,a);
36: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
37: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
39: MatCreate(MPI_COMM_SELF,size,size,size,size,&B);
40: MatSetType(B,MATSEQDENSE);
41: MatSeqDenseSetPreallocation(B,b);
42: MatSeqDenseSetLDA(B,lda);
43: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
44: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
46: PetscMalloc(size*sizeof(PetscScalar),&x);
47: for (i=0; i<size; i++) {
48: x[i] = one;
49: }
50: VecCreateSeqWithArray(MPI_COMM_SELF,size,x,&X);
51: VecAssemblyBegin(X);
52: VecAssemblyEnd(X);
54: PetscMalloc(size*sizeof(PetscScalar),&y);
55: VecCreateSeqWithArray(MPI_COMM_SELF,size,y,&Y);
56: VecAssemblyBegin(Y);
57: VecAssemblyEnd(Y);
59: PetscMalloc(size*sizeof(PetscScalar),&z);
60: VecCreateSeqWithArray(MPI_COMM_SELF,size,z,&Z);
61: VecAssemblyBegin(Z);
62: VecAssemblyEnd(Z);
64: /*
65: * Solve with A and B
66: */
67: KSPCreate(MPI_COMM_SELF,&solver);
68: KSPSetType(solver,KSPPREONLY);
69: KSPGetPC(solver,&pc);
70: PCSetType(pc,PCLU);
71: KSPSetOperators(solver,A,A,DIFFERENT_NONZERO_PATTERN);
72: KSPSetRhs(solver,X);
73: KSPSetSolution(solver,Y);
74: KSPSolve(solver);
75: KSPSetOperators(solver,B,B,DIFFERENT_NONZERO_PATTERN);
76: KSPSetSolution(solver,Z);
77: KSPSolve(solver);
78: VecAXPY(&mone,Y,Z);
79: VecNorm(Z,NORM_2,&nrm);
80: printf("Test1; error norm=%e\n",nrm);
82: PetscFinalize();
83: return 0;
84: }