1
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";
9
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
22
23 rhs
24 ---
25 0 3.0
26 1 12.0
27 2 6.0
28
29 solu
30 ----
31 0 1.0
32 0 1.0
33 0 1.0
34 */
35
36
37 #include <petscmat.h>
38
main(int argc,char ** args)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];
44 PetscErrorCode ierr;
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;
53
54 ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
55 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
56 if (size != 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"This is a uniprocessor example only!");
57
58 /* Read in matrix, rhs and exact solution from ascii files */
59 ierr = PetscOptionsGetString(NULL,NULL,"-Ain",Ain,sizeof(Ain),&flg_A);CHKERRQ(ierr);
60 ierr = PetscOptionsHasName(NULL,NULL,"-noshift",&flg);CHKERRQ(ierr);
61 if (flg) shift = 0;
62 if (flg_A) {
63 ierr = PetscPrintf(PETSC_COMM_SELF,"\n Read matrix in ascii format ...\n");CHKERRQ(ierr);
64 ierr = PetscFOpen(PETSC_COMM_SELF,Ain,"r",&Afile);CHKERRQ(ierr);
65 nsizes = 3;
66 ierr = PetscOptionsGetIntArray(NULL,NULL,"-nosizesinfile",sizes,&nsizes,&flg);CHKERRQ(ierr);
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 ierr = PetscPrintf(PETSC_COMM_SELF,"m: %d, n: %d, nz: %d \n", m,n,nz);CHKERRQ(ierr);
76 if (m != n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ, "Number of rows, cols must be same for this example\n");
77 ierr = MatCreate(PETSC_COMM_SELF,&A);CHKERRQ(ierr);
78 ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m,n);CHKERRQ(ierr);
79 ierr = MatSetFromOptions(A);CHKERRQ(ierr);
80 ierr = MatSeqAIJSetPreallocation(A,nz/m,NULL);CHKERRQ(ierr);
81 ierr = MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);CHKERRQ(ierr);
82
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 ierr = MatSetValues(A,1,&row,1,&col,&val,INSERT_VALUES);CHKERRQ(ierr);
87 }
88 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
89 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
90 fclose(Afile);
91 }
92
93 ierr = PetscOptionsGetString(NULL,NULL,"-rhs",rhs,sizeof(rhs),&flg_b);CHKERRQ(ierr);
94 if (flg_b) {
95 ierr = VecCreate(PETSC_COMM_SELF,&b);CHKERRQ(ierr);
96 ierr = VecSetSizes(b,PETSC_DECIDE,n);CHKERRQ(ierr);
97 ierr = VecSetFromOptions(b);CHKERRQ(ierr);
98 ierr = PetscPrintf(PETSC_COMM_SELF,"\n Read rhs in ascii format ...\n");CHKERRQ(ierr);
99 ierr = PetscFOpen(PETSC_COMM_SELF,rhs,"r",&bfile);CHKERRQ(ierr);
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 ierr = VecSetValues(b,1,&i,&val,INSERT_VALUES);CHKERRQ(ierr);
103 }
104 ierr = VecAssemblyBegin(b);CHKERRQ(ierr);
105 ierr = VecAssemblyEnd(b);CHKERRQ(ierr);
106 fclose(bfile);
107 }
108
109 ierr = PetscOptionsGetString(NULL,NULL,"-solu",solu,sizeof(solu),&flg_u);CHKERRQ(ierr);
110 if (flg_u) {
111 ierr = VecCreate(PETSC_COMM_SELF,&u);CHKERRQ(ierr);
112 ierr = VecSetSizes(u,PETSC_DECIDE,n);CHKERRQ(ierr);
113 ierr = VecSetFromOptions(u);CHKERRQ(ierr);
114 ierr = PetscPrintf(PETSC_COMM_SELF,"\n Read exact solution in ascii format ...\n");CHKERRQ(ierr);
115 ierr = PetscFOpen(PETSC_COMM_SELF,solu,"r",&ufile);CHKERRQ(ierr);
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 ierr = VecSetValues(u,1,&i,&val,INSERT_VALUES);CHKERRQ(ierr);
119 }
120 ierr = VecAssemblyBegin(u);CHKERRQ(ierr);
121 ierr = VecAssemblyEnd(u);CHKERRQ(ierr);
122 fclose(ufile);
123 }
124
125 /* Write matrix, rhs and exact solution in Petsc binary file */
126 ierr = PetscPrintf(PETSC_COMM_SELF,"\n Write matrix in binary to 'matrix.dat' ...\n");CHKERRQ(ierr);
127 ierr = PetscViewerBinaryOpen(PETSC_COMM_SELF,"matrix.dat",FILE_MODE_WRITE,&view);CHKERRQ(ierr);
128 ierr = MatView(A,view);CHKERRQ(ierr);
129
130 if (flg_b) { /* Write rhs in Petsc binary file */
131 ierr = PetscPrintf(PETSC_COMM_SELF,"\n Write rhs in binary to 'matrix.dat' ...\n");CHKERRQ(ierr);
132 ierr = VecView(b,view);CHKERRQ(ierr);
133 }
134 if (flg_u) {
135 ierr = PetscPrintf(PETSC_COMM_SELF,"\n Write exact solution in binary to 'matrix.dat' ...\n");CHKERRQ(ierr);
136 ierr = VecView(u,view);CHKERRQ(ierr);
137 }
138
139 /* Check accuracy of the data */
140 if (flg_A & flg_b & flg_u) {
141 ierr = VecDuplicate(u,&u_tmp);CHKERRQ(ierr);
142 ierr = MatMult(A,u,u_tmp);CHKERRQ(ierr);
143 ierr = VecAXPY(u_tmp,-1.0,b);CHKERRQ(ierr);
144 ierr = VecNorm(u_tmp,NORM_2,&res_norm);CHKERRQ(ierr);
145 ierr = PetscPrintf(PETSC_COMM_SELF,"\n Accuracy of the reading data: | b - A*u |_2 : %g \n",res_norm);CHKERRQ(ierr);
146 ierr = VecDestroy(&u_tmp);CHKERRQ(ierr);
147 }
148
149 ierr = MatDestroy(&A);CHKERRQ(ierr);
150 if (flg_b) {ierr = VecDestroy(&b);CHKERRQ(ierr);}
151 if (flg_u) {ierr = VecDestroy(&u);CHKERRQ(ierr);}
152 ierr = PetscViewerDestroy(&view);CHKERRQ(ierr);
153 ierr = PetscFinalize();
154 return ierr;
155 }
156
157 /*TEST
158
159 build:
160 requires: !define(PETSC_USE_64BIT_INDICES) double !complex
161
162 test:
163 requires: datafilespath
164 args: -Ain ${DATAFILESPATH}/matrices/indefinite/afiro_A.dat
165
166 TEST*/
167