1 /*  testExtract.c  */
2 
3 #include "../InpMtx.h"
4 #include "../../Drand.h"
5 
6 /*--------------------------------------------------------------------*/
7 int
main(int argc,char * argv[])8 main ( int argc, char *argv[] )
9 /*
10    ------------------------------------------------------------------
11    generate a random matrix and test the submatrix extraction method.
12    the output is a matlab file to test correctness.
13 
14    created -- 98oct15, cca
15  --------------------------------------------------------------------
16 */
17 {
18 int        coordType, dataType, ii, msglvl, ncolA, ncol1, ncol2,
19            nitem, nrowA, nrow1, nrow2, rc, seed, symmetryflag ;
20 int        *colids, *rowids, *temp ;
21 InpMtx     *A ;
22 IV         *cols1IV, *cols2IV, *rows1IV, *rows2IV ;
23 FILE       *msgFile ;
24 
25 if ( argc != 14 ) {
26    fprintf(stdout,
27    "\n\n %% usage : %s msglvl msgFile dataType symmetryflag coordType "
28    "\n %%         nrowA ncolA nitem nrow1 nrow2 ncol1 ncol2 seed"
29    "\n %%    msglvl   -- message level"
30    "\n %%    msgFile  -- message file"
31    "\n %%    dataType -- type of matrix entries"
32    "\n %%       1 -- real"
33    "\n %%       2 -- complex"
34    "\n %%    symmetryflag -- type of matrix symmetry"
35    "\n %%       0 -- symmetric"
36    "\n %%       1 -- hermitian"
37    "\n %%       2 -- nonsymmetric"
38    "\n %%    coordType -- coordinate Type"
39    "\n %%       1 -- by rows"
40    "\n %%       2 -- by columns"
41    "\n %%       3 -- by chevrons"
42    "\n %%    nrowA    -- number of rows in A"
43    "\n %%    ncolA    -- number of columns in A"
44    "\n %%    nitem    -- number of items to be loaded into A"
45    "\n %%    nrow1    -- number of rows in the first block"
46    "\n %%    nrow2    -- number of rows in the second block"
47    "\n %%    ncol1    -- number of columns in the first block"
48    "\n %%    ncol2    -- number of columns in the second block"
49    "\n %%    seed     -- random number seed"
50    "\n", argv[0]) ;
51    return(0) ;
52 }
53 msglvl = atoi(argv[1]) ;
54 if ( strcmp(argv[2], "stdout") == 0 ) {
55    msgFile = stdout ;
56 } else if ( (msgFile = fopen(argv[2], "a")) == NULL ) {
57    fprintf(stderr, "\n fatal error in %s"
58            "\n unable to open file %s\n",
59            argv[0], argv[2]) ;
60    return(-1) ;
61 }
62 dataType      = atoi(argv[3]) ;
63 symmetryflag  = atoi(argv[4]) ;
64 coordType     = atoi(argv[5]) ;
65 nrowA         = atoi(argv[6]) ;
66 ncolA         = atoi(argv[7]) ;
67 nitem         = atoi(argv[8]) ;
68 nrow1         = atoi(argv[9]) ;
69 nrow2         = atoi(argv[10]) ;
70 ncol1         = atoi(argv[11]) ;
71 ncol2         = atoi(argv[12]) ;
72 seed          = atoi(argv[13]) ;
73 fprintf(msgFile,
74         "\n %% %s "
75         "\n %% msglvl        -- %d"
76         "\n %% msgFile       -- %s"
77         "\n %% dataType      -- %d"
78         "\n %% symmetryflag  -- %d"
79         "\n %% coordType     -- %d"
80         "\n %% nrowA         -- %d"
81         "\n %% ncolA         -- %d"
82         "\n %% nitem         -- %d"
83         "\n %% nrow1         -- %d"
84         "\n %% nrow2         -- %d"
85         "\n %% ncol1         -- %d"
86         "\n %% ncol2         -- %d"
87         "\n %% seed          -- %d"
88         "\n",
89         argv[0], msglvl, argv[2], dataType, symmetryflag, coordType,
90         nrowA, ncolA, nitem, nrow1, nrow2, ncol1, ncol2, seed) ;
91 fflush(msgFile) ;
92 if ( dataType != SPOOLES_REAL && dataType != SPOOLES_COMPLEX ) {
93    fprintf(stderr, "\n invalid value %d for dataType\n", dataType) ;
94    exit(-1) ;
95 }
96 if (  symmetryflag != SPOOLES_SYMMETRIC
97    && symmetryflag != SPOOLES_HERMITIAN
98    && symmetryflag != SPOOLES_NONSYMMETRIC ) {
99    fprintf(stderr,
100            "\n invalid value %d for symmetryflag\n", symmetryflag) ;
101    exit(-1) ;
102 }
103 if ( symmetryflag == SPOOLES_HERMITIAN && dataType != SPOOLES_COMPLEX ){
104    fprintf(stderr,
105            "\n symmetryflag is hermitian, data type is not complex\n") ;
106    exit(-1) ;
107 }
108 if (  coordType != INPMTX_BY_ROWS
109    && coordType != INPMTX_BY_COLUMNS
110    && coordType != INPMTX_BY_CHEVRONS ) {
111    fprintf(stderr,
112            "\n invalid value %d for coordType\n", coordType) ;
113    exit(-1) ;
114 }
115 if ( nrowA <= 0 || ncolA <= 0 || nitem <= 0 ) {
116    fprintf(stderr,
117            "\n invalid value: nrow = %d, ncol = %d, nitem = %d",
118            nrowA, ncolA, nitem) ;
119    exit(-1) ;
120 }
121 fprintf(msgFile, "\n %% symmetryflag %d, nrowA %d, ncolA %d\n",
122         symmetryflag, nrowA, ncolA) ;
123 if (   (   symmetryflag == SPOOLES_SYMMETRIC
124         || symmetryflag == SPOOLES_HERMITIAN)
125     && nrowA != ncolA ) {
126    fprintf(stderr,
127            "\n symmetric matrix, nrowA %d, ncolA %d\n", nrowA, ncolA) ;
128    exit(-1) ;
129 }
130 if ( nrow1 < 0 || nrow2 < 0 || (nrow1 + nrow2 != nrowA) ) {
131    fprintf(stderr,
132            "\n invalid value: nrow = %d, nrow1 = %d, nrow2 = %d",
133            nrowA, nrow1, nrow2) ;
134    exit(-1) ;
135 }
136 if ( ncol1 < 0 || ncol2 < 0 || (ncol1 + ncol2 != ncolA) ) {
137    fprintf(stderr,
138            "\n invalid value: ncol = %d, ncol1 = %d, ncol2 = %d",
139            ncolA, ncol1, ncol2) ;
140    exit(-1) ;
141 }
142 /*
143    ----------------------------
144    initialize the matrix object
145    ----------------------------
146 */
147 A = InpMtx_new() ;
148 rc = InpMtx_randomMatrix(A, dataType, coordType, INPMTX_BY_VECTORS,
149                    nrowA, ncolA, symmetryflag, 0, nitem, seed) ;
150 /*
151    -------------------------------------------
152    write the assembled matrix to a matlab file
153    -------------------------------------------
154 */
155 fprintf(msgFile, "\n A = zeros(%d,%d) ;", nrowA, ncolA) ;
156 InpMtx_writeForMatlab(A, "A", msgFile) ;
157 if ( symmetryflag == SPOOLES_SYMMETRIC ) {
158    fprintf(msgFile,
159           "\n A = diag(diag(A)) + triu(A,1) + transpose(triu(A,1)) ;") ;
160 } else if ( symmetryflag == SPOOLES_HERMITIAN ) {
161    fprintf(msgFile,
162          "\n A = diag(diag(A)) + triu(A,1) + ctranspose(triu(A,1)) ;") ;
163 }
164 /*
165    ----------------------------------
166    generate row and column id vectors
167    ----------------------------------
168 */
169 temp = IVinit(nrowA, -1) ;
170 IVramp(nrowA, temp, 0, 1) ;
171 IVshuffle(nrowA, temp, ++seed) ;
172 if ( nrow1 > 0 ) {
173    rows1IV = IV_new() ;
174    IV_init(rows1IV, nrow1, NULL) ;
175    IVqsortUp(nrow1, temp) ;
176    IVcopy(nrow1, IV_entries(rows1IV), temp) ;
177 } else {
178    rows1IV = NULL ;
179 }
180 if ( nrow2 > 0 ) {
181    rows2IV = IV_new() ;
182    IV_init(rows2IV, nrow2, NULL) ;
183    IVqsortUp(nrow2, temp + nrow1) ;
184    IVcopy(nrow2, IV_entries(rows2IV), temp + nrow1) ;
185 } else {
186    rows2IV = NULL ;
187 }
188 IVfree(temp) ;
189 if ( symmetryflag == SPOOLES_NONSYMMETRIC ) {
190    temp = IVinit(ncolA, -1) ;
191    IVramp(ncolA, temp, 0, 1) ;
192    IVshuffle(ncolA, temp, ++seed) ;
193    if ( ncol1 > 0 ) {
194       cols1IV = IV_new() ;
195       IV_init(cols1IV, ncol1, NULL) ;
196       IVqsortUp(ncol1, temp) ;
197       IVcopy(ncol1, IV_entries(cols1IV), temp) ;
198    } else {
199       cols1IV = NULL ;
200    }
201    if ( ncol2 > 0 ) {
202       cols2IV = IV_new() ;
203       IV_init(cols2IV, ncol2, NULL) ;
204       IVqsortUp(ncol2, temp + ncol1) ;
205       IVcopy(ncol2, IV_entries(cols2IV), temp + ncol1) ;
206    } else {
207       cols2IV = NULL ;
208    }
209    IVfree(temp) ;
210 } else {
211    if ( ncol1 > 0 ) {
212       cols1IV = IV_new() ;
213       IV_init(cols1IV, ncol1, NULL) ;
214       IV_copy(cols1IV, rows1IV) ;
215    } else {
216       cols1IV = NULL ;
217    }
218    if ( ncol2 > 0 ) {
219       cols2IV = IV_new() ;
220       IV_init(cols2IV, ncol2, NULL) ;
221       IV_copy(cols2IV, rows2IV) ;
222    } else {
223       cols2IV = NULL ;
224    }
225 }
226 if ( nrow1 > 0 ) {
227    fprintf(msgFile, "\n rows1 = zeros(%d,1) ;", nrow1) ;
228    IV_writeForMatlab(rows1IV, "rows1", msgFile) ;
229 }
230 if ( nrow2 > 0 ) {
231    fprintf(msgFile, "\n rows2 = zeros(%d,1) ;", nrow2) ;
232    IV_writeForMatlab(rows2IV, "rows2", msgFile) ;
233 }
234 if ( ncol1 > 0 ) {
235    fprintf(msgFile, "\n cols1 = zeros(%d,1) ;", ncol1) ;
236    IV_writeForMatlab(cols1IV, "cols1", msgFile) ;
237 }
238 if ( ncol2 > 0 ) {
239    fprintf(msgFile, "\n cols2 = zeros(%d,1) ;", ncol2) ;
240    IV_writeForMatlab(cols2IV, "cols2", msgFile) ;
241 }
242 /*
243    -----------------------
244    extract the submatrices
245    -----------------------
246 */
247 if ( nrow1 > 0 ) {
248    if ( ncol1 > 0 ) {
249       InpMtx   *A11 = InpMtx_new() ;
250       rc = InpMtx_initFromSubmatrix(A11, A, rows1IV, cols1IV,
251                                    symmetryflag, msglvl, msgFile) ;
252       if ( rc != 1 ) {
253          fprintf(stderr, "\n error return %d for A11\n", rc) ;
254          exit(-1) ;
255       }
256       fprintf(msgFile, "\n A11 = zeros(%d,%d) ;", nrow1, ncol1) ;
257       InpMtx_writeForMatlab(A11, "A11", msgFile) ;
258       if ( symmetryflag == SPOOLES_SYMMETRIC ) {
259          fprintf(msgFile,
260   "\n A11 = diag(diag(A11)) + triu(A11,1) + transpose(triu(A11,1)) ;") ;
261       } else if ( symmetryflag == SPOOLES_HERMITIAN ) {
262          fprintf(msgFile,
263  "\n A11 = diag(diag(A11)) + triu(A11,1) + ctranspose(triu(A11,1)) ;") ;
264       }
265       fprintf(msgFile,
266               "\n err11 = max(max(abs(A11 - A(rows1,cols1)))) ;") ;
267       InpMtx_free(A11) ;
268    } else {
269       fprintf(msgFile, "\n err11 = 0 ;") ;
270    }
271    if ( ncol2 > 0 ) {
272       InpMtx   *A12 = InpMtx_new() ;
273       rc = InpMtx_initFromSubmatrix(A12, A, rows1IV, cols2IV,
274                                    symmetryflag, msglvl, msgFile) ;
275       if ( rc != 1 ) {
276          fprintf(stderr, "\n error return %d for A12\n", rc) ;
277          exit(-1) ;
278       }
279       fprintf(msgFile, "\n A12 = zeros(%d,%d) ;", nrow1, ncol2) ;
280       InpMtx_writeForMatlab(A12, "A12", msgFile) ;
281       fprintf(msgFile,
282               "\n err12 = max(max(abs(A12 - A(rows1,cols2)))) ;") ;
283       InpMtx_free(A12) ;
284    } else {
285       fprintf(msgFile, "\n err12 = 0 ;") ;
286    }
287 } else {
288    fprintf(msgFile, "\n err11 = 0 ;") ;
289    fprintf(msgFile, "\n err12 = 0 ;") ;
290 }
291 if ( nrow2 > 0 ) {
292    if ( ncol1 > 0 ) {
293       InpMtx   *A21 = InpMtx_new() ;
294       rc = InpMtx_initFromSubmatrix(A21, A, rows2IV, cols1IV,
295                                    symmetryflag, msglvl, msgFile) ;
296       if ( rc != 1 ) {
297          fprintf(stderr, "\n error return %d for A21\n", rc) ;
298          exit(-1) ;
299       }
300       fprintf(msgFile, "\n A21 = zeros(%d,%d) ;", nrow2, ncol1) ;
301       InpMtx_writeForMatlab(A21, "A21", msgFile) ;
302       fprintf(msgFile,
303               "\n err21 = max(max(abs(A21 - A(rows2,cols1)))) ;") ;
304       InpMtx_free(A21) ;
305    } else {
306       fprintf(msgFile, "\n err21 = 0 ;") ;
307    }
308    if ( ncol2 > 0 ) {
309       InpMtx   *A22 = InpMtx_new() ;
310       rc = InpMtx_initFromSubmatrix(A22, A, rows2IV, cols2IV,
311                                    symmetryflag, msglvl, msgFile) ;
312       if ( rc != 1 ) {
313          fprintf(stderr, "\n error return %d for A22\n", rc) ;
314          exit(-1) ;
315       }
316       fprintf(msgFile, "\n A22 = zeros(%d,%d) ;", nrow2, ncol2) ;
317       InpMtx_writeForMatlab(A22, "A22", msgFile) ;
318       if ( symmetryflag == SPOOLES_SYMMETRIC ) {
319          fprintf(msgFile,
320   "\n A22 = diag(diag(A22)) + triu(A22,1) + transpose(triu(A22,1)) ;") ;
321       } else if ( symmetryflag == SPOOLES_HERMITIAN ) {
322          fprintf(msgFile,
323  "\n A22 = diag(diag(A22)) + triu(A22,1) + ctranspose(triu(A22,1)) ;") ;
324       }
325       fprintf(msgFile,
326               "\n err22 = max(max(abs(A22 - A(rows2,cols2)))) ;") ;
327       InpMtx_free(A22) ;
328    } else {
329       fprintf(msgFile, "\n err22 = 0 ;") ;
330    }
331 } else {
332    fprintf(msgFile, "\n err21 = 0 ;") ;
333    fprintf(msgFile, "\n err22 = 0 ;") ;
334 }
335 fprintf(msgFile, "\n error = [ err11 err12 ; err21 err22 ] ") ;
336 /*
337    ------------------------
338    free the working storage
339    ------------------------
340 */
341 InpMtx_free(A) ;
342 if ( rows1IV != NULL ) { IV_free(rows1IV) ; }
343 if ( rows2IV != NULL ) { IV_free(rows2IV) ; }
344 if ( cols1IV != NULL ) { IV_free(cols1IV) ; }
345 if ( cols2IV != NULL ) { IV_free(cols2IV) ; }
346 
347 fclose(msgFile) ;
348 
349 return(1) ; }
350 
351 /*--------------------------------------------------------------------*/
352