Actual source code: ex78.c
petsc-3.14.6 2021-03-30
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: -rhs <filename> : input rhs in ascii format\n\
8: -solu <filename> : input true solution in ascii format\n\\n";
10: /*
11: Example: ./ex78 -Ain Ain -rhs rhs -solu solu -noshift -mat_view
12: with the datafiles in the followig format:
13: Ain (I and J start at 0):
14: ------------------------
15: 3 3 6
16: 0 0 1.0
17: 0 1 2.0
18: 1 0 3.0
19: 1 1 4.0
20: 1 2 5.0
21: 2 2 6.0
23: rhs
24: ---
25: 0 3.0
26: 1 12.0
27: 2 6.0
29: solu
30: ----
31: 0 1.0
32: 0 1.0
33: 0 1.0
34: */
37: #include <petscmat.h>
39: int main(int argc,char **args)
40: {
41: Mat A = NULL;
42: Vec b,u = NULL,u_tmp;
43: char Ain[PETSC_MAX_PATH_LEN],rhs[PETSC_MAX_PATH_LEN],solu[PETSC_MAX_PATH_LEN];
45: int m,n = 0,nz,dummy; /* these are fscaned so kept as int */
46: PetscInt i,col,row,shift = 1,sizes[3],nsizes;
47: PetscScalar val;
48: PetscReal res_norm;
49: FILE *Afile,*bfile,*ufile;
50: PetscViewer view;
51: PetscBool flg_A,flg_b,flg_u,flg;
52: PetscMPIInt size;
54: PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
55: MPI_Comm_size(PETSC_COMM_WORLD,&size);
56: if (size != 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"This is a uniprocessor example only!");
58: /* Read in matrix, rhs and exact solution from ascii files */
59: PetscOptionsGetString(NULL,NULL,"-Ain",Ain,sizeof(Ain),&flg_A);
60: PetscOptionsHasName(NULL,NULL,"-noshift",&flg);
61: if (flg) shift = 0;
62: if (flg_A) {
63: PetscPrintf(PETSC_COMM_SELF,"\n Read matrix in ascii format ...\n");
64: PetscFOpen(PETSC_COMM_SELF,Ain,"r",&Afile);
65: nsizes = 3;
66: PetscOptionsGetIntArray(NULL,NULL,"-nosizesinfile",sizes,&nsizes,&flg);
67: if (flg) {
68: if (nsizes != 3) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"Must pass in three m,n,nz as arguments for -nosizesinfile");
69: m = sizes[0];
70: n = sizes[1];
71: nz = sizes[2];
72: } else {
73: if (fscanf(Afile,"%d %d %d\n",&m,&n,&nz) != 3) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Badly formatted input file\n");
74: }
75: PetscPrintf(PETSC_COMM_SELF,"m: %d, n: %d, nz: %d \n", m,n,nz);
76: if (m != n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ, "Number of rows, cols must be same for this example\n");
77: MatCreate(PETSC_COMM_SELF,&A);
78: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m,n);
79: MatSetFromOptions(A);
80: MatSeqAIJSetPreallocation(A,nz/m,NULL);
81: MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
83: for (i=0; i<nz; i++) {
84: if (fscanf(Afile,"%d %d %le\n",&row,&col,(double*)&val) != 3) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Badly formatted input file\n");
85: row -= shift; col -= shift; /* set index set starts at 0 */
86: MatSetValues(A,1,&row,1,&col,&val,INSERT_VALUES);
87: }
88: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
89: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
90: fclose(Afile);
91: }
93: PetscOptionsGetString(NULL,NULL,"-rhs",rhs,sizeof(rhs),&flg_b);
94: if (flg_b) {
95: VecCreate(PETSC_COMM_SELF,&b);
96: VecSetSizes(b,PETSC_DECIDE,n);
97: VecSetFromOptions(b);
98: PetscPrintf(PETSC_COMM_SELF,"\n Read rhs in ascii format ...\n");
99: PetscFOpen(PETSC_COMM_SELF,rhs,"r",&bfile);
100: for (i=0; i<n; i++) {
101: if (fscanf(bfile,"%d %le\n",&dummy,(double*)&val) != 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Badly formatted input file\n");
102: VecSetValues(b,1,&i,&val,INSERT_VALUES);
103: }
104: VecAssemblyBegin(b);
105: VecAssemblyEnd(b);
106: fclose(bfile);
107: }
109: PetscOptionsGetString(NULL,NULL,"-solu",solu,sizeof(solu),&flg_u);
110: if (flg_u) {
111: VecCreate(PETSC_COMM_SELF,&u);
112: VecSetSizes(u,PETSC_DECIDE,n);
113: VecSetFromOptions(u);
114: PetscPrintf(PETSC_COMM_SELF,"\n Read exact solution in ascii format ...\n");
115: PetscFOpen(PETSC_COMM_SELF,solu,"r",&ufile);
116: for (i=0; i<n; i++) {
117: if (fscanf(ufile,"%d %le\n",&dummy,(double*)&val) != 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Badly formatted input file\n");
118: VecSetValues(u,1,&i,&val,INSERT_VALUES);
119: }
120: VecAssemblyBegin(u);
121: VecAssemblyEnd(u);
122: fclose(ufile);
123: }
125: /* Write matrix, rhs and exact solution in Petsc binary file */
126: PetscPrintf(PETSC_COMM_SELF,"\n Write matrix in binary to 'matrix.dat' ...\n");
127: PetscViewerBinaryOpen(PETSC_COMM_SELF,"matrix.dat",FILE_MODE_WRITE,&view);
128: MatView(A,view);
130: if (flg_b) { /* Write rhs in Petsc binary file */
131: PetscPrintf(PETSC_COMM_SELF,"\n Write rhs in binary to 'matrix.dat' ...\n");
132: VecView(b,view);
133: }
134: if (flg_u) {
135: PetscPrintf(PETSC_COMM_SELF,"\n Write exact solution in binary to 'matrix.dat' ...\n");
136: VecView(u,view);
137: }
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: PetscPrintf(PETSC_COMM_SELF,"\n Accuracy of the reading data: | b - A*u |_2 : %g \n",res_norm);
146: VecDestroy(&u_tmp);
147: }
149: MatDestroy(&A);
150: if (flg_b) {VecDestroy(&b);}
151: if (flg_u) {VecDestroy(&u);}
152: PetscViewerDestroy(&view);
153: PetscFinalize();
154: return ierr;
155: }
157: /*TEST
159: build:
160: requires: !define(PETSC_USE_64BIT_INDICES) double !complex
162: test:
163: requires: datafilespath
164: args: -Ain ${DATAFILESPATH}/matrices/indefinite/afiro_A.dat
166: TEST*/