1 /* testInput.c */
2
3 #include "../../timings.h"
4 #include "../InpMtx.h"
5
6 /*--------------------------------------------------------------------*/
7 int
main(int argc,char * argv[])8 main ( int argc, char *argv[] )
9 /*
10 -----------------------------------------------
11 test the InpMtx input methods for Peter Schartz
12
13 created -- 98sep04, cca
14 -----------------------------------------------
15 */
16 {
17 double estpar, growth, t1, t2 ;
18 double *entries ;
19 FILE *msgFile ;
20 int count, ii, irow, maxsize, msglvl, nent, neqns, n1, n2, n3,
21 size, size1, size2, type ;
22 int *indices, *indices1, *indices2, *list ;
23 InpMtx *mtxA ;
24 IVL *adjIVL, *fullIVL, *lowerIVL ;
25
26 if ( argc != 9 ) {
27 fprintf(stdout,
28 "\n\n usage : %s msglvl msgFile n1 n2 n3 estpar growth"
29 "\n msglvl -- message level"
30 "\n msgFile -- message file"
31 "\n type -- type of entries"
32 "\n 0 -- indices only"
33 "\n 1 -- real entries"
34 "\n 2 -- complex entries"
35 "\n n1 -- # of grid points in first direction"
36 "\n n2 -- # of grid points in second direction"
37 "\n n3 -- # of grid points in third direction"
38 "\n estpar -- estimation for nent"
39 "\n growth -- growth factor"
40 "\n", argv[0]) ;
41 return(0) ;
42 }
43 msglvl = atoi(argv[1]) ;
44 if ( strcmp(argv[2], "stdout") == 0 ) {
45 msgFile = stdout ;
46 } else if ( (msgFile = fopen(argv[2], "a")) == NULL ) {
47 fprintf(stderr, "\n fatal error in %s"
48 "\n unable to open file %s\n",
49 argv[0], argv[2]) ;
50 return(-1) ;
51 }
52 type = atoi(argv[3]) ;
53 n1 = atoi(argv[4]) ;
54 n2 = atoi(argv[5]) ;
55 n3 = atoi(argv[6]) ;
56 estpar = atof(argv[7]) ;
57 growth = atof(argv[8]) ;
58 fprintf(msgFile,
59 "\n %s "
60 "\n msglvl -- %d"
61 "\n msgFile -- %s"
62 "\n type -- %d"
63 "\n n1 -- %d"
64 "\n n2 -- %d"
65 "\n n3 -- %d"
66 "\n estpar -- %f"
67 "\n growth -- %f"
68 "\n",
69 argv[0], msglvl, argv[2], type, n1, n2, n3, estpar, growth) ;
70 fflush(msgFile) ;
71 if ( n1 <= 0 || n2 <= 0 || n3 <= 0 || estpar < 0.0 || growth <= 1.0 ) {
72 fprintf(stderr, "\n fatal error in testInput, bad input\n") ;
73 exit(-1) ;
74 }
75 /*
76 -----------------------------------
77 set up the grid adjacency structure
78 -----------------------------------
79 */
80 neqns = n1 * n2 * n3 ;
81 MARKTIME(t1) ;
82 if ( n1 == 1 ) {
83 adjIVL = IVL_make9P(n2, n3, 1) ;
84 } else if ( n2 == 1 ) {
85 adjIVL = IVL_make9P(n1, n3, 1) ;
86 } else if ( n3 == 1 ) {
87 adjIVL = IVL_make9P(n1, n2, 1) ;
88 } else {
89 adjIVL = IVL_make27P(n1, n2, n3, 1) ;
90 }
91 nent = IVL_tsize(adjIVL) ;
92 MARKTIME(t2) ;
93 fprintf(msgFile, "\n\n CPU %8.3f : make full adjacency, %d entries",
94 t2 - t1, nent) ;
95 if ( msglvl > 1 ) {
96 fprintf(msgFile, "\n\n full adjacency structure, %d entries", nent) ;
97 IVL_writeForHumanEye(adjIVL, msgFile) ;
98 }
99 /*
100 ----------------------------------
101 make the lower adjacency structure
102 ----------------------------------
103 */
104 MARKTIME(t1) ;
105 lowerIVL = IVL_new() ;
106 IVL_init1(lowerIVL, IVL_CHUNKED, neqns) ;
107 list = IVinit(neqns, -1) ;
108 for ( irow = 0 ; irow < neqns ; irow++ ) {
109 IVL_listAndSize(adjIVL, irow, &size, &indices) ;
110 for ( ii = count = 0 ; ii < size ; ii++ ) {
111 if ( indices[ii] >= irow ) {
112 list[count++] = indices[ii] ;
113 }
114 }
115 IVL_setList(lowerIVL, irow, count, list) ;
116 }
117 nent = IVL_tsize(lowerIVL) ;
118 MARKTIME(t2) ;
119 fprintf(msgFile, "\n\n CPU %8.3f : make lower adjacency, %d entries",
120 t2 - t1, nent) ;
121 if ( msglvl > 1 ) {
122 fprintf(msgFile, "\n\n lower adjacency structure, %d entries", nent);
123 IVL_writeForHumanEye(adjIVL, msgFile) ;
124 }
125 /*
126 ---------------------------------------------------
127 create a vector to hold entries,
128 its size is the maximum size of the lower adjacency
129 ---------------------------------------------------
130 */
131 maxsize = IVL_maxListSize(adjIVL) ;
132 entries = DVinit(2*maxsize, 0.0) ;
133 /*
134 ----------------------------
135 initialize the InpMtx object
136 ----------------------------
137 */
138 MARKTIME(t1) ;
139 mtxA = InpMtx_new() ;
140 InpMtx_init(mtxA, INPMTX_BY_COLUMNS, type, estpar*nent, 0) ;
141 InpMtx_setResizeMultiple(mtxA, growth) ;
142 MARKTIME(t2) ;
143 fprintf(msgFile, "\n\n CPU %8.3f : initialize InpMtx", t2 - t1) ;
144
145 /*
146 ----------------
147 load the columns
148 ----------------
149 */
150 MARKTIME(t1) ;
151 if ( INPMTX_IS_INDICES_ONLY(mtxA) ) {
152 for ( irow = 0 ; irow < neqns ; irow++ ) {
153 IVL_listAndSize(lowerIVL, irow, &size, &indices) ;
154 InpMtx_inputColumn(mtxA, irow, size, indices) ;
155 }
156 } else if ( INPMTX_IS_REAL_ENTRIES(mtxA) ) {
157 for ( irow = 0 ; irow < neqns ; irow++ ) {
158 IVL_listAndSize(lowerIVL, irow, &size, &indices) ;
159 InpMtx_inputRealColumn(mtxA, irow, size, indices, entries) ;
160 }
161 } else if ( INPMTX_IS_COMPLEX_ENTRIES(mtxA) ) {
162 for ( irow = 0 ; irow < neqns ; irow++ ) {
163 IVL_listAndSize(lowerIVL, irow, &size, &indices) ;
164 InpMtx_inputComplexColumn(mtxA, irow, size, indices, entries) ;
165 }
166 }
167 MARKTIME(t2) ;
168 fprintf(msgFile, "\n\n CPU %8.3f : load entries by columns", t2 - t1) ;
169 if ( msglvl > 1 ) {
170 fprintf(msgFile, "\n\n mtxA") ;
171 InpMtx_writeForHumanEye(mtxA, msgFile) ;
172 }
173 /*
174 -----------------------------
175 sort and compress the entries
176 -----------------------------
177 */
178 MARKTIME(t1) ;
179 InpMtx_sortAndCompress(mtxA) ;
180 MARKTIME(t2) ;
181 fprintf(msgFile, "\n\n CPU %8.3f : sort and compress", t2 - t1) ;
182 /*
183 -------------------
184 set the vector mode
185 -------------------
186 */
187 MARKTIME(t1) ;
188 InpMtx_changeStorageMode(mtxA, INPMTX_BY_VECTORS) ;
189 MARKTIME(t2) ;
190 fprintf(msgFile, "\n\n CPU %8.3f : convert to vectors", t2 - t1) ;
191 /*
192 --------------------------------------
193 construct the full adjacency structure
194 --------------------------------------
195 */
196 MARKTIME(t1) ;
197 fullIVL = InpMtx_fullAdjacency(mtxA) ;
198 MARKTIME(t2) ;
199 fprintf(msgFile, "\n\n CPU %8.3f : construct the full adjacency",
200 t2 - t1) ;
201 /*
202 -----------------------------------------
203 compare the two full adjacency structures
204 -----------------------------------------
205 */
206 for ( irow = 0 ; irow < neqns ; irow++ ) {
207 IVL_listAndSize(adjIVL, irow, &size1, &indices1) ;
208 IVL_listAndSize(fullIVL, irow, &size2, &indices2) ;
209 if ( size1 != size2 ) {
210 fprintf(msgFile, "\n\n error, irow %d, size1 %d, size2 %d",
211 irow, size1, size2) ;
212 exit(-1) ;
213 }
214 for ( ii = 0 ; ii < size1 ; ii++ ) {
215 if ( indices1[ii] != indices2[ii] ) {
216 fprintf(msgFile, "\n\n error, irow %d", irow) ;
217 fprintf(msgFile, "\n indices1") ;
218 IVfprintf(msgFile, size1, indices1) ;
219 fprintf(msgFile, "\n indices2") ;
220 IVfprintf(msgFile, size1, indices2) ;
221 exit(-1) ;
222 }
223 }
224 }
225 /*
226 ------------------------
227 free the working storage
228 ------------------------
229 */
230 IVL_free(adjIVL) ;
231 IVL_free(lowerIVL) ;
232 IVL_free(fullIVL) ;
233 InpMtx_free(mtxA) ;
234 DVfree(entries) ;
235 IVfree(list) ;
236
237 fprintf(msgFile, "\n") ;
238 fclose(msgFile) ;
239
240 return(1) ; }
241
242 /*--------------------------------------------------------------------*/
243