Actual source code: ex78.c
2: 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\
3: Writes them using the PETSc sparse format.\n\
4: Note: I and J start at 1, not 0, use -noshift if indices in file start with zero!\n\
5: Input parameters are:\n\
6: -Ain <filename> : input matrix in ascii format\n\
7: -bin <filename> : input rhs in ascii format\n\
8: -uin <filename> : input true solution in ascii format\n\
9: Run this program: ex33h -Ain Ain -bin bin -uin uin\n\n";
11: #include petscmat.h
15: int main(int argc,char **args)
16: {
17: Mat A;
18: Vec b,u,u_tmp;
19: char Ain[PETSC_MAX_PATH_LEN],bin[PETSC_MAX_PATH_LEN],uin[PETSC_MAX_PATH_LEN];
21: int m,n,nz,dummy,*col=0,*row=0; /* these are fscaned so kept as int */
22: PetscInt i,col_i,row_i,*nnz,*ib,shift = 1,sizes[3],nsizes;
23: PetscScalar val_i,*work=0;
24: PetscReal res_norm,*val=0,*bval=0,*uval=0;
25: FILE *Afile,*bfile,*ufile;
26: PetscViewer view;
27: PetscTruth flg_A,flg_b,flg_u,flg;
29: PetscInitialize(&argc,&args,(char *)0,help);
31: /* Read in matrix, rhs and exact solution from an ascii file */
32: PetscOptionsGetString(PETSC_NULL,"-Ain",Ain,PETSC_MAX_PATH_LEN-1,&flg_A);
33: PetscOptionsHasName(PETSC_NULL,"-noshift",&flg);
34: if (flg) shift = 0;
35: if (flg_A){
36: PetscPrintf(PETSC_COMM_SELF,"\n Read matrix in ascii format ...\n");
37: PetscFOpen(PETSC_COMM_SELF,Ain,"r",&Afile);
38: nsizes = 3;
39: PetscOptionsGetIntArray(PETSC_NULL,"-nosizesinfile",sizes,&nsizes,&flg);
40: if (flg) {
41: if (nsizes != 3) SETERRQ(1,"Must pass in three m,n,nz as arguments for -nosizesinfile");
42: m = sizes[0];
43: n = sizes[1];
44: nz = sizes[2];
45: } else {
46: fscanf(Afile,"%d %d %d\n",&m,&n,&nz);
47: }
48: printf("m: %d, n: %d, nz: %d \n", m,n,nz);
49: if (m != n) SETERRQ(PETSC_ERR_ARG_SIZ, "Number of rows, cols must be same for SBAIJ format\n");
51: VecCreateSeq(PETSC_COMM_SELF,n,&b);
52: VecCreateSeq(PETSC_COMM_SELF,n,&u);
54: PetscMalloc(m*sizeof(PetscInt),&nnz);
55: PetscMemzero(nnz,m*sizeof(PetscInt));
56: PetscMalloc(nz*sizeof(PetscReal),&val);
57: PetscMalloc(nz*sizeof(PetscInt),&row);
58: PetscMalloc(nz*sizeof(PetscInt),&col);
59: for (i=0; i<nz; i++) {
60: fscanf(Afile,"%d %d %le\n",row+i,col+i,(double*)val+i); /* modify according to data file! */
61: row[i] -= shift; col[i] -= shift; /* set index set starts at 0 */
62: nnz[row[i]]++;
63: }
64: printf("row[0]: %d, col[0]: %d, val: %g\n",row[0],col[0],val[0]);
65: printf("row[last]: %d, col: %d, val: %g\n",row[nz-1],col[nz-1],val[nz-1]);
66: fflush(stdout);
67: fclose(Afile);
68: MatCreateSeqBAIJ(PETSC_COMM_SELF,1,n,n,0,nnz,&A);
69: PetscFree(nnz);
70: }
72: PetscOptionsGetString(PETSC_NULL,"-bin",bin,PETSC_MAX_PATH_LEN-1,&flg_b);
73: if (flg_b){
74: PetscPrintf(PETSC_COMM_SELF,"\n Read rhs in ascii format ...\n");
75: PetscFOpen(PETSC_COMM_SELF,bin,"r",&bfile);
76: PetscMalloc(n*sizeof(PetscReal),&bval);
77: PetscMalloc(n*sizeof(PetscScalar),&work);
78: PetscMalloc(n*sizeof(PetscInt),&ib);
79: for (i=0; i<n; i++) {
80: /* fscanf(bfile,"%d %le\n",ib+i,bval+i); ib[i]--; */ /* modify according to data file! */
81: fscanf(bfile,"%d %le\n",&dummy,(double*)(bval+i)); ib[i] = i; /* modify according to data file! */
82: }
83: printf("bval[0]: %g, bval[%d]: %g\n",bval[0],(int)ib[n-1],bval[n-1]);
84: fflush(stdout);
85: fclose(bfile);
86: }
88: PetscOptionsGetString(PETSC_NULL,"-uin",uin,PETSC_MAX_PATH_LEN-1,&flg_u);
89: if (flg_u){
90: PetscPrintf(PETSC_COMM_SELF,"\n Read exact solution in ascii format ...\n");
91: PetscFOpen(PETSC_COMM_SELF,uin,"r",&ufile);
92: PetscMalloc(n*sizeof(PetscReal),&uval);
93: for (i=0; i<n; i++) {
94: fscanf(ufile,"%d %le\n",&dummy,(double*)(uval+i)); /* modify according to data file! */
95: }
96: printf("uval[0]: %g, uval[%d]: %g\n",uval[0], n-1, uval[n-1]);
97: fflush(stdout);
98: fclose(ufile);
99: }
101: if(flg_A){
102: /*
103: for (i=0; i<n; i++){
104: MatSetValues(A,1,&i,1,&i,&zero,INSERT_VALUES);
105: }
106: */
107: for (i=0; i<nz; i++) {
108: row_i =(int)row[i]; col_i =(int)col[i]; val_i = (PetscScalar)val[i];
109: MatSetValues(A,1,&row_i,1,&col_i,&val_i,INSERT_VALUES);
110: }
111: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
112: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
113: PetscFree(col);
114: PetscFree(val);
115: PetscFree(row);
116: }
117: if(flg_b){
118: for (i=0; i<n; i++) work[i]=(PetscScalar)bval[i];
119: VecSetValues(b,n,ib,work,INSERT_VALUES);
120: VecAssemblyBegin(b);
121: VecAssemblyEnd(b);
122: /* printf("b: \n"); VecView(b,PETSC_VIEWER_STDOUT_SELF); */
123: PetscFree(bval);
124: }
126: if(flg_u & flg_b){
127: for (i=0; i<n; i++) work[i]=(PetscScalar)uval[i];
128: VecSetValues(u,n,ib,work,INSERT_VALUES);
129: VecAssemblyBegin(u);
130: VecAssemblyEnd(u);
131: /* printf("u: \n"); VecView(u,PETSC_VIEWER_STDOUT_SELF); */
132: PetscFree(uval);
133: }
134:
135: if(flg_b) {
136: PetscFree(ib);
137: PetscFree(work);
138: }
139: /* Check accuracy of the data */
140: if (flg_A & flg_b & flg_u){
141: VecDuplicate(u,&u_tmp);
142: MatMult(A,u,u_tmp);
143: VecAXPY(u_tmp,-1.0,b);
144: VecNorm(u_tmp,NORM_2,&res_norm);
145: printf("\n Accuracy of the reading data: | b - A*u |_2 : %g \n",res_norm);
147: /* Write the matrix, rhs and exact solution in Petsc binary file */
148: PetscPrintf(PETSC_COMM_SELF,"\n Write matrix and rhs in binary to 'matrix.dat' ...\n");
149: PetscViewerBinaryOpen(PETSC_COMM_SELF,"matrix.dat",FILE_MODE_WRITE,&view);
150: MatView(A,view);
151: VecView(b,view);
152: VecView(u,view);
153: PetscViewerDestroy(view);
155: VecDestroy(b);
156: VecDestroy(u);
157: VecDestroy(u_tmp);
158:
159: MatDestroy(A);
160: }
162: PetscFinalize();
163: return 0;
164: }