1 /* test_maxabs.c */
2
3 #include "../Chv.h"
4 #include "../../Drand.h"
5 #include "../../timings.h"
6
7 /*--------------------------------------------------------------------*/
8 int
main(int argc,char * argv[])9 main ( int argc, char *argv[] )
10 /*
11 ----------------------------------------------------
12 test the Chv_maxabsInDiagonal(), Chv_maxabsInRow()
13 and Chv_maxabsInColumn() methods.
14 the program's output is a matlab file
15 to check correctness of the code.
16
17 created -- 98apr22, cca
18 ----------------------------------------------------
19 */
20 {
21 Chv *chv ;
22 double imag, maxval, real, t1, t2 ;
23 double *colmaxes, *entries, *rowmaxes ;
24 Drand *drand ;
25 FILE *msgFile ;
26 int icol, ii, irow, jcol, jrow, msglvl, ncol, nD, nent,
27 nL, nrow, nU, rc, seed, symflag, tag, type ;
28 int *colind, *colmark, *rowind, *rowmark ;
29
30 if ( argc != 8 ) {
31 fprintf(stdout,
32 "\n\n usage : %s msglvl msgFile nD nU type symflag seed "
33 "\n msglvl -- message level"
34 "\n msgFile -- message file"
35 "\n nD -- # of rows and columns in the (1,1) block"
36 "\n nU -- # of columns in the (1,2) block"
37 "\n type -- entries type"
38 "\n 1 --> real"
39 "\n 2 --> complex"
40 "\n symflag -- symmetry flag"
41 "\n 0 --> symmetric"
42 "\n 1 --> hermitian"
43 "\n 2 --> nonsymmetric "
44 "\n seed -- random number seed"
45 "\n", argv[0]) ;
46 return(0) ;
47 }
48 if ( (msglvl = atoi(argv[1])) < 0 ) {
49 fprintf(stderr, "\n message level must be positive\n") ;
50 exit(-1) ;
51 }
52 if ( strcmp(argv[2], "stdout") == 0 ) {
53 msgFile = stdout ;
54 } else if ( (msgFile = fopen(argv[2], "a")) == NULL ) {
55 fprintf(stderr, "\n unable to open file %s\n", argv[2]) ;
56 return(-1) ;
57 }
58 nD = atoi(argv[3]) ;
59 nU = atoi(argv[4]) ;
60 type = atoi(argv[5]) ;
61 symflag = atoi(argv[6]) ;
62 seed = atoi(argv[7]) ;
63 fprintf(msgFile, "\n %% testChv:"
64 "\n %% msglvl = %d"
65 "\n %% msgFile = %s"
66 "\n %% nD = %d"
67 "\n %% nU = %d"
68 "\n %% type = %d"
69 "\n %% symflag = %d"
70 "\n %% seed = %d",
71 msglvl, argv[2], nD, nU, type, symflag, seed) ;
72 nL = nU ;
73 nrow = nD + nL ;
74 ncol = nD + nU ;
75 /*
76 -----------------------------
77 check for errors in the input
78 -----------------------------
79 */
80 if ( nD <= 0 || nL < 0 || nU < 0
81 || symflag < 0 || symflag > 3 ) {
82 fprintf(stderr, "\n invalid input"
83 "\n nD = %d, nL = %d, nU = %d, symflag = %d\n",
84 nD, nL, nU, symflag) ;
85 exit(-1) ;
86 }
87 fprintf(msgFile,
88 "\n nD = %d ;"
89 "\n nL = %d ;"
90 "\n nU = %d ;"
91 "\n nrow = nD + nL ;"
92 "\n ncol = nD + nU ;",
93 nD, nL, nU) ;
94 /*
95 --------------------------------------
96 initialize the random number generator
97 --------------------------------------
98 */
99 drand = Drand_new() ;
100 Drand_init(drand) ;
101 Drand_setSeed(drand, seed) ;
102 Drand_setNormal(drand, 0.0, 1.0) ;
103 /*
104 ----------------------------
105 initialize the Chv object
106 ----------------------------
107 */
108 MARKTIME(t1) ;
109 chv = Chv_new() ;
110 Chv_init(chv, 0, nD, nL, nU, type, symflag) ;
111 MARKTIME(t2) ;
112 fprintf(msgFile, "\n %% CPU : %.3f to initialize chv object",
113 t2 - t1) ;
114 fflush(msgFile) ;
115 Chv_columnIndices(chv, &ncol, &colind) ;
116 IVramp(ncol, colind, 0, 1) ;
117 if ( CHV_IS_NONSYMMETRIC(chv) ) {
118 Chv_rowIndices(chv, &nrow, &rowind) ;
119 IVramp(nrow, rowind, 0, 1) ;
120 }
121 /*
122 ------------------------------------
123 load the entries with random entries
124 ------------------------------------
125 */
126 nent = Chv_nent(chv) ;
127 entries = Chv_entries(chv) ;
128 if ( CHV_IS_REAL(chv) ) {
129 Drand_fillDvector(drand, nent, entries) ;
130 } else if ( CHV_IS_COMPLEX(chv) ) {
131 Drand_fillDvector(drand, 2*nent, entries) ;
132 }
133 if ( CHV_IS_HERMITIAN(chv) ) {
134 /*
135 ---------------------------------------------------------
136 hermitian example, set imaginary part of diagonal to zero
137 ---------------------------------------------------------
138 */
139 for ( irow = 0 ; irow < nD ; irow++ ) {
140 Chv_complexEntry(chv, irow, irow, &real, &imag) ;
141 Chv_setComplexEntry(chv, irow, irow, real, 0.0) ;
142 }
143 }
144 fprintf(msgFile, "\n %% matrix entries") ;
145 Chv_writeForMatlab(chv, "a", msgFile) ;
146 /*
147 -----------------------------
148 find the row and column maxes
149 -----------------------------
150 */
151 rowmaxes = DVinit(nrow, 0.0) ;
152 colmaxes = DVinit(nrow, 0.0) ;
153 for ( irow = 0 ; irow < nD ; irow++ ) {
154 jcol = Chv_maxabsInRow(chv, irow, &maxval) ;
155 rowmaxes[irow] = maxval ;
156 }
157 for ( jcol = 0 ; jcol < nD ; jcol++ ) {
158 irow = Chv_maxabsInColumn(chv, jcol, &maxval) ;
159 colmaxes[jcol] = maxval ;
160 }
161 fprintf(msgFile, "\n\n rowmaxes = [ ...") ;
162 for ( irow = 0 ; irow < nD ; irow++ ) {
163 fprintf(msgFile, "\n %20.12e", rowmaxes[irow]) ;
164 }
165 fprintf(msgFile, " ] ;") ;
166 fprintf(msgFile, "\n\n colmaxes = [ ...") ;
167 for ( jcol = 0 ; jcol < nD ; jcol++ ) {
168 fprintf(msgFile, "\n %20.12e", colmaxes[jcol]) ;
169 }
170 fprintf(msgFile, " ] ;") ;
171 /*
172 -----------------
173 check with matlab
174 -----------------
175 */
176 fprintf(msgFile,
177 "\n\n for irow = 1:nD"
178 "\n maxval = max(abs(a(irow,:))) ;"
179 "\n rmaxes(irow) = maxval ;"
180 "\nend"
181 "\nrowerror = norm(rmaxes' - rowmaxes) "
182 "\nfor jcol = 1:nD"
183 "\n maxval = max(abs(a(:,jcol))) ;"
184 "\n cmaxes(jcol) = maxval ;"
185 "\nend"
186 "\ncolerror = norm(cmaxes' - colmaxes) ") ;
187 /*
188 -----------------------------------------------------
189 find the row and column maxes of just the (1,1) block
190 -----------------------------------------------------
191 */
192 rowmark = IVinit(nD, -1) ;
193 colmark = IVinit(nD, -1) ;
194 tag = -1 ;
195 for ( irow = 0 ; irow < nD ; irow++ ) {
196 jcol = Chv_maxabsInRow11(chv, irow, colmark, tag, &maxval) ;
197 rowmaxes[irow] = maxval ;
198 }
199 for ( jcol = 0 ; jcol < nD ; jcol++ ) {
200 irow = Chv_maxabsInColumn11(chv, jcol, rowmark, tag, &maxval) ;
201 colmaxes[jcol] = maxval ;
202 }
203 fprintf(msgFile, "\n\n rowmaxes = [ ...") ;
204 for ( irow = 0 ; irow < nD ; irow++ ) {
205 fprintf(msgFile, "\n %20.12e", rowmaxes[irow]) ;
206 }
207 fprintf(msgFile, " ] ;") ;
208 fprintf(msgFile, "\n\n colmaxes = [ ...") ;
209 for ( jcol = 0 ; jcol < nD ; jcol++ ) {
210 fprintf(msgFile, "\n %20.12e", colmaxes[jcol]) ;
211 }
212 fprintf(msgFile, " ] ;") ;
213 /*
214 -----------------
215 check with matlab
216 -----------------
217 */
218 fprintf(msgFile,
219 "\n\n for irow = 1:nD"
220 "\n maxval = max(abs(a(irow,1:nD))) ;"
221 "\n rmaxes(irow) = maxval ;"
222 "\nend"
223 "\nrow11error = norm(rmaxes' - rowmaxes) "
224 "\nfor jcol = 1:nD"
225 "\n maxval = max(abs(a(1:nD,jcol))) ;"
226 "\n cmaxes(jcol) = maxval ;"
227 "\nend"
228 "\ncol11error = norm(cmaxes' - colmaxes) ") ;
229 /*
230 ---------------------------------------------
231 find the diagonal max of just the (1,1) block
232 ---------------------------------------------
233 */
234 jcol = Chv_maxabsInDiagonal11(chv, colmark, tag, &maxval) ;
235 fprintf(msgFile, "\n\n maxval = %20.12e ;", maxval) ;
236 /*
237 -----------------
238 check with matlab
239 -----------------
240 */
241 fprintf(msgFile,
242 "\n\n maxabs = abs(a(1,1)) ;"
243 "\nfor irow = 2:nD"
244 "\n val = abs(a(irow,irow)) ;"
245 "\n if maxabs < val"
246 "\n maxabs = val ;"
247 "\n end "
248 "\nend"
249 "\ndiag11error = abs(maxabs - maxval)") ;
250 fprintf(msgFile,
251 "\n [ rowerror colerror row11error col11error diag11error]") ;
252 /*
253 ------------------------
254 free the working storage
255 ------------------------
256 */
257 Chv_free(chv) ;
258 Drand_free(drand) ;
259 DVfree(rowmaxes) ;
260 DVfree(colmaxes) ;
261 IVfree(rowmark) ;
262 IVfree(colmark) ;
263
264 fprintf(msgFile, "\n") ;
265
266 return(1) ; }
267
268 /*--------------------------------------------------------------------*/
269