Actual source code: ex78.c
1: /*$Id: ex78.c,v 1.14 2001/08/07 21:30:08 bsmith Exp $*/
3: static char help[] = "Reads in a matrix in ASCII Matlab format (I,J,A), read in vectors rhs and exact_solu in ASCII format.\n\
4: Writes them using the PETSc sparse format.\n\
5: Note: I and J start at 1, not 0!\n\
6: Input parameters are:\n\
7: -Ain <filename> : input matrix in ascii format\n\
8: -bin <filename> : input rhs in ascii format\n\
9: -uin <filename> : input true solution in ascii format\n\
10: Run this program: ex33h -Ain Ain -bin bin -uin uin\n\n";
12: #include petscmat.h
16: int main(int argc,char **args)
17: {
18: Mat A;
19: Vec b,u,u_tmp;
20: char Ain[128],bin[128],uin[128];
21: int i,m,n,nz,ierr,*ib=0,col_i,row_i;
22: PetscScalar val_i,*work=0,mone=-1.0;
23: PetscReal *col=0,*row=0,res_norm,*val=0,*bval=0,*uval=0;
24: FILE *Afile,*bfile,*ufile;
25: PetscViewer view;
26: PetscTruth flg_A,flg_b,flg_u;
28: PetscInitialize(&argc,&args,(char *)0,help);
30: /* Read in matrix, rhs and exact solution from an ascii file */
31: PetscOptionsGetString(PETSC_NULL,"-Ain",Ain,127,&flg_A);
32: if (flg_A){
33: PetscPrintf(PETSC_COMM_SELF,"\n Read matrix in ascii format ...\n");
34: PetscFOpen(PETSC_COMM_SELF,Ain,"r",&Afile);
35: fscanf(Afile,"%d %d %d\n",&m,&n,&nz);
36: printf("m: %d, n: %d, nz: %d \n", m,n,nz);
37: if (m != n) SETERRQ(PETSC_ERR_ARG_SIZ, "Number of rows, cols must be same for SBAIJ format\n");
39: MatCreateSeqBAIJ(PETSC_COMM_SELF,1,n,n,20,0,&A);
40: VecCreateSeq(PETSC_COMM_SELF,n,&b);
41: VecCreateSeq(PETSC_COMM_SELF,n,&u);
43: PetscMalloc(nz*sizeof(PetscReal),&val);
44: PetscMalloc(nz*sizeof(PetscReal),&row);
45: PetscMalloc(nz*sizeof(PetscReal),&col);
46: for (i=0; i<nz; i++) {
47: fscanf(Afile,"%le %le %le\n",row+i,col+i,val+i); /* modify according to data file! */
48: row[i]--; col[i]--; /* set index set starts at 0 */
49: }
50: printf("row[0]: %g, col[0]: %g, val: %g\n",row[0],col[0],val[0]);
51: printf("row[last]: %g, col: %g, val: %g\n",row[nz-1],col[nz-1],val[nz-1]);
52: fclose(Afile);
53: }
55: PetscOptionsGetString(PETSC_NULL,"-bin",bin,127,&flg_b);
56: if (flg_b){
57: PetscPrintf(PETSC_COMM_SELF,"\n Read rhs in ascii format ...\n");
58: PetscFOpen(PETSC_COMM_SELF,bin,"r",&bfile);
59: PetscMalloc(n*sizeof(PetscReal),&bval);
60: PetscMalloc(n*sizeof(PetscScalar),&work);
61: PetscMalloc(n*sizeof(int),&ib);
62: for (i=0; i<n; i++) {
63: /* fscanf(bfile,"%d %le\n",ib+i,bval+i); ib[i]--; */ /* modify according to data file! */
64: fscanf(bfile,"%le\n",bval+i); ib[i] = i; /* modify according to data file! */
65: }
66: printf("bval[0]: %g, bval[%d]: %g\n",bval[0],ib[n-1],bval[n-1]);
67: fclose(bfile);
68: }
70: PetscOptionsGetString(PETSC_NULL,"-uin",uin,127,&flg_u);
71: if (flg_u){
72: PetscPrintf(PETSC_COMM_SELF,"\n Read exact solution in ascii format ...\n");
73: PetscFOpen(PETSC_COMM_SELF,uin,"r",&ufile);
74: PetscMalloc(n*sizeof(PetscReal),&uval);
75: for (i=0; i<n; i++) {
76: fscanf(ufile," %le\n",uval+i); /* modify according to data file! */
77: }
78: printf("uval[0]: %g, uval[%d]: %g\n",uval[0], n-1, uval[n-1]);
79: fclose(ufile);
80: }
82: if(flg_A){
83: /*
84: for (i=0; i<n; i++){
85: MatSetValues(A,1,&i,1,&i,&zero,INSERT_VALUES);
86: }
87: */
88: for (i=0; i<nz; i++) {
89: row_i =(int)row[i]; col_i =(int)col[i]; val_i = (PetscScalar)val[i];
90: MatSetValues(A,1,&row_i,1,&col_i,&val_i,INSERT_VALUES);
91: }
92: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
93: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
94: PetscFree(col);
95: PetscFree(val);
96: PetscFree(row);
97: }
98: if(flg_b){
99: for (i=0; i<n; i++) work[i]=(PetscScalar)bval[i];
100: VecSetValues(b,n,ib,work,INSERT_VALUES);
101: VecAssemblyBegin(b);
102: VecAssemblyEnd(b);
103: /* printf("b: \n"); VecView(b,PETSC_VIEWER_STDOUT_SELF); */
104: PetscFree(bval);
105: }
107: if(flg_u & flg_b){
108: for (i=0; i<n; i++) work[i]=(PetscScalar)uval[i];
109: VecSetValues(u,n,ib,work,INSERT_VALUES);
110: VecAssemblyBegin(u);
111: VecAssemblyEnd(u);
112: /* printf("u: \n"); VecView(u,PETSC_VIEWER_STDOUT_SELF); */
113: PetscFree(uval);
114: }
115:
116: if(flg_b) {
117: PetscFree(ib);
118: PetscFree(work);
119: }
120: /* Check accuracy of the data */
121: if (flg_A & flg_b & flg_u){
122: VecDuplicate(u,&u_tmp);
123: MatMult(A,u,u_tmp);
124: VecAXPY(&mone,b,u_tmp);
125: VecNorm(u_tmp,NORM_2,&res_norm);
126: printf("\n Accuracy of the reading data: | b - A*u |_2 : %g \n",res_norm);
128: /* Write the matrix, rhs and exact solution in Petsc binary file */
129: PetscPrintf(PETSC_COMM_SELF,"\n Write matrix and rhs in binary to 'matrix.dat' ...\n");
130: PetscViewerBinaryOpen(PETSC_COMM_SELF,"matrix.dat",PETSC_BINARY_CREATE,&view);
131: MatView(A,view);
132: VecView(b,view);
133: VecView(u,view);
134: PetscViewerDestroy(view);
136: VecDestroy(b);
137: VecDestroy(u);
138: VecDestroy(u_tmp);
139:
140: MatDestroy(A);
141: }
143: PetscFinalize();
144: return 0;
145: }