1 /* extract.c */
2
3 #include "../InpMtx.h"
4
5 /*--------------------------------------------------------------------*/
6 /*
7 ----------------------------------------------------------
8 purpose -- to extract a submatrix, B = A(BrowsIV, BcolsIV)
9
10 return values ---
11 1 -- normal return
12 -1 -- B is NULL
13 -2 -- BcolsIV is NULL
14 -3 -- BrowsIV is NULL
15 -4 -- B is NULL
16 -5 -- invalid input mode for A
17 -6 -- invalid coordinate type for A
18 -7 -- invalid symmetryflag
19 -8 -- hermitian flag but not complex
20 -9 -- msglvl > 0 and msgFile = NULL
21
22 created -- 98oct15, cca
23 ----------------------------------------------------------
24 */
25 int
InpMtx_initFromSubmatrix(InpMtx * B,InpMtx * A,IV * BrowsIV,IV * BcolsIV,int symmetryflag,int msglvl,FILE * msgFile)26 InpMtx_initFromSubmatrix (
27 InpMtx *B,
28 InpMtx *A,
29 IV *BrowsIV,
30 IV *BcolsIV,
31 int symmetryflag,
32 int msglvl,
33 FILE *msgFile
34 ) {
35 double *dbuf, *entA ;
36 int colA, colB, ii, jj, jjfirst, jjlast, kk, maxcolA, maxn,
37 maxrowA, nbuf, ncolB, nent, nrowB, nvector, rowA, rowB,
38 oldCoordType, oldStorageMode, rowsAndColumnsAreIdentical ;
39 int *colmap, *colsA, *colsB, *ibuf1, *ibuf2, *offsets, *rowmap,
40 *rowsB, *sizes, *vecids ;
41 /*
42 ---------------
43 check the input
44 ---------------
45 */
46 if ( B == NULL ) {
47 fprintf(stderr, "\n error in InpMtx_initFromSubmatrix()"
48 "\n B is NULL\n") ;
49 return(-1) ;
50 }
51 if ( BrowsIV == NULL ) {
52 fprintf(stderr, "\n error in InpMtx_initFromSubmatrix()"
53 "\n BrowsIV is NULL\n") ;
54 return(-2) ;
55 }
56 if ( BcolsIV == NULL ) {
57 fprintf(stderr, "\n error in InpMtx_initFromSubmatrix()"
58 "\n BcolsIV is NULL\n") ;
59 return(-3) ;
60 }
61 if ( A == NULL ) {
62 fprintf(stderr, "\n error in InpMtx_initFromSubmatrix()"
63 "\n A is NULL\n") ;
64 return(-4) ;
65 }
66 switch ( A->inputMode ) {
67 case SPOOLES_REAL :
68 case SPOOLES_COMPLEX :
69 case INPMTX_INDICES_ONLY :
70 break ;
71 default :
72 fprintf(stderr, "\n error in InpMtx_initFromSubmatrix()"
73 "\n invalid inputMode %d for A\n", A->inputMode) ;
74 return(-5) ;
75 break ;
76 }
77 switch ( A->coordType ) {
78 case INPMTX_BY_ROWS :
79 case INPMTX_BY_COLUMNS :
80 case INPMTX_BY_CHEVRONS :
81 break ;
82 default :
83 fprintf(stderr, "\n error in InpMtx_initFromSubmatrix()"
84 "\n invalid coordType %d for A\n", A->coordType) ;
85 return(-6) ;
86 break ;
87 }
88 switch ( symmetryflag ) {
89 case SPOOLES_SYMMETRIC :
90 case SPOOLES_HERMITIAN :
91 case SPOOLES_NONSYMMETRIC :
92 break ;
93 default :
94 fprintf(stderr, "\n error in InpMtx_initFromSubmatrix()"
95 "\n invalid symmetryflag %d \n", symmetryflag) ;
96 return(-7) ;
97 break ;
98 }
99 if ( A->inputMode != SPOOLES_COMPLEX
100 && symmetryflag == SPOOLES_HERMITIAN ) {
101 fprintf(stderr, "\n error in InpMtx_initFromSubmatrix()"
102 "\n Hermitian symmetry flag but A is not complex\n") ;
103 return(-8) ;
104 }
105 if ( msglvl > 0 && msgFile == NULL ) {
106 fprintf(stderr, "\n error in InpMtx_initFromSubmatrix()"
107 "\n msglvl = %d, msgFile = NULL\n", msglvl) ;
108 return(-9) ;
109 }
110 /*--------------------------------------------------------------------*/
111 /*
112 ------------
113 initialize B
114 ------------
115 */
116 InpMtx_init(B, INPMTX_BY_ROWS, A->inputMode, 0, 0) ;
117 /*
118 -----------------------------
119 get the range of entries in A
120 -----------------------------
121 */
122 InpMtx_range(A, NULL, &maxcolA, NULL, &maxrowA) ;
123 maxn = (maxcolA >= maxrowA) ? maxcolA : maxrowA ;
124 /*
125 --------------------------------------------------
126 get the # and indices of the rows and columns of B
127 --------------------------------------------------
128 */
129 IV_sizeAndEntries(BrowsIV, &nrowB, &rowsB) ;
130 IV_sizeAndEntries(BcolsIV, &ncolB, &colsB) ;
131 if ( nrowB != ncolB ) {
132 rowsAndColumnsAreIdentical = 0 ;
133 } else {
134 rowsAndColumnsAreIdentical = 1 ;
135 for ( ii = 0 ; ii < nrowB ; ii++ ) {
136 if ( rowsB[ii] != colsB[ii] ) {
137 rowsAndColumnsAreIdentical = 0 ;
138 break ;
139 }
140 }
141 }
142 /*
143 ---------------------------
144 set up the local column map
145 ---------------------------
146 */
147 colmap = IVinit(1+maxn, -1) ;
148 for ( ii = 0 ; ii < ncolB ; ii++ ) {
149 colmap[colsB[ii]] = ii ;
150 }
151 if ( msglvl > 2 ) {
152 fprintf(msgFile, "\n colmap") ;
153 IVfprintf(msgFile, 1+maxn, colmap) ;
154 fflush(msgFile) ;
155 }
156 /*
157 ---------------------------------------
158 change the coordinate type of A to rows
159 and storage mode to vectors
160 ---------------------------------------
161 */
162 if ( (oldCoordType = A->coordType) != INPMTX_BY_ROWS ) {
163 InpMtx_changeCoordType(A, INPMTX_BY_ROWS) ;
164 }
165 if ( (oldStorageMode = A->storageMode) != INPMTX_BY_VECTORS ) {
166 InpMtx_changeStorageMode(A, INPMTX_BY_VECTORS) ;
167 }
168 if ( msglvl > 2 ) {
169 fprintf(msgFile, "\n A") ;
170 InpMtx_writeForHumanEye(A, msgFile) ;
171 fflush(msgFile) ;
172 }
173 /*
174 -----------------------------------------------------------
175 switch over the input mode. search the rows of A that will
176 belong to B, load into the buffer all entries that belong
177 in columns of B. when the buffer is full, add to B matrix.
178 -----------------------------------------------------------
179 */
180 nbuf = 100 ;
181 ibuf1 = IVinit(nbuf, -1) ;
182 ibuf2 = IVinit(nbuf, -1) ;
183 kk = 0 ;
184 if ( INPMTX_IS_REAL_ENTRIES(A) ) {
185 dbuf = DVinit(nbuf, 0.0) ;
186 for ( rowB = 0 ; rowB < nrowB ; rowB++ ) {
187 rowA = rowsB[rowB] ;
188 InpMtx_realVector(A, rowA, &nent, &colsA, &entA) ;
189 for ( jj = 0 ; jj < nent ; jj++ ) {
190 if ( (colB = colmap[colsA[jj]]) != -1 ) {
191 ibuf1[kk] = rowB ; ibuf2[kk] = colB ;
192 dbuf[kk] = entA[jj] ;
193 if ( ++kk == nbuf ) {
194 InpMtx_inputRealTriples(B, kk, ibuf1, ibuf2, dbuf) ;
195 kk = 0 ;
196 }
197 }
198 }
199 }
200 } else if ( INPMTX_IS_COMPLEX_ENTRIES(A) ) {
201 dbuf = DVinit(2*nbuf, 0.0) ;
202 for ( rowB = 0 ; rowB < nrowB ; rowB++ ) {
203 rowA = rowsB[rowB] ;
204 InpMtx_complexVector(A, rowA, &nent, &colsA, &entA) ;
205 for ( jj = 0 ; jj < nent ; jj++ ) {
206 if ( (colB = colmap[colsA[jj]]) != -1 ) {
207 ibuf1[kk] = rowB ; ibuf2[kk] = colB ;
208 dbuf[2*kk] = entA[2*jj] ; dbuf[2*kk+1] = entA[2*jj+1] ;
209 if ( ++kk == nbuf ) {
210 InpMtx_inputComplexTriples(B, kk, ibuf1, ibuf2, dbuf) ;
211 kk = 0 ;
212 }
213 }
214 }
215 }
216 } else if ( INPMTX_IS_INDICES_ONLY(A) ) {
217 dbuf = NULL ;
218 for ( rowB = 0 ; rowB < nrowB ; rowB++ ) {
219 rowA = rowsB[rowB] ;
220 InpMtx_vector(A, rowA, &nent, &colsA) ;
221 for ( jj = 0 ; jj < nent ; jj++ ) {
222 if ( (colB = colmap[colsA[jj]]) != -1 ) {
223 ibuf1[kk] = rowB ; ibuf2[kk] = colB ;
224 if ( ++kk == nbuf ) {
225 InpMtx_inputTriples(B, kk, ibuf1, ibuf2) ;
226 kk = 0 ;
227 }
228 }
229 }
230 }
231 }
232 if ( ( symmetryflag == SPOOLES_SYMMETRIC
233 || symmetryflag == SPOOLES_HERMITIAN)
234 && !rowsAndColumnsAreIdentical ) {
235 /*
236 ----------------------------------------------
237 matrix is symmetric or hermitian, search lower
238 triangle of A for entries that belong to B.
239 ----------------------------------------------
240 */
241 rowmap = IVinit(1+maxn, -1) ;
242 for ( ii = 0 ; ii < nrowB ; ii++ ) {
243 rowmap[rowsB[ii]] = ii ;
244 }
245 if ( msglvl > 2 ) {
246 fprintf(msgFile, "\n rowmap") ;
247 IVfprintf(msgFile, 1+maxn, rowmap) ;
248 fflush(msgFile) ;
249 }
250 if ( INPMTX_IS_REAL_ENTRIES(A) ) {
251 for ( colB = 0 ; colB < ncolB ; colB++ ) {
252 colA = colsB[colB] ;
253 InpMtx_realVector(A, colA, &nent, &colsA, &entA) ;
254 for ( jj = 0 ; jj < nent ; jj++ ) {
255 if ( (rowB = rowmap[colsA[jj]]) != -1 ) {
256 ibuf1[kk] = rowB ; ibuf2[kk] = colB ;
257 dbuf[kk] = entA[jj] ;
258 if ( ++kk == nbuf ) {
259 InpMtx_inputRealTriples(B, kk, ibuf1, ibuf2, dbuf) ;
260 kk = 0 ;
261 }
262 }
263 }
264 }
265 } else if ( INPMTX_IS_COMPLEX_ENTRIES(A) ) {
266 for ( colB = 0 ; colB < ncolB ; colB++ ) {
267 colA = colsB[colB] ;
268 InpMtx_complexVector(A, colA, &nent, &colsA, &entA) ;
269 for ( jj = 0 ; jj < nent ; jj++ ) {
270 if ( (rowB = rowmap[colsA[jj]]) != -1 ) {
271 ibuf1[kk] = rowB ; ibuf2[kk] = colB ;
272 dbuf[2*kk] = entA[2*jj] ; dbuf[2*kk+1] = -entA[2*jj+1] ;
273 if ( ++kk == nbuf ) {
274 InpMtx_inputComplexTriples(B, kk, ibuf1, ibuf2, dbuf);
275 kk = 0 ;
276 }
277 }
278 }
279 }
280 } else if ( INPMTX_IS_INDICES_ONLY(A) ) {
281 for ( colB = 0 ; colB < ncolB ; colB++ ) {
282 colA = colsB[colB] ;
283 InpMtx_vector(A, colA, &nent, &colsA) ;
284 for ( jj = 0 ; jj < nent ; jj++ ) {
285 if ( (rowB = rowmap[colsA[jj]]) != -1 ) {
286 ibuf1[kk] = rowB ; ibuf2[kk] = colB ;
287 if ( ++kk == nbuf ) {
288 InpMtx_inputTriples(B, kk, ibuf1, ibuf2) ;
289 kk = 0 ;
290 }
291 }
292 }
293 }
294 }
295 IVfree(rowmap) ;
296 }
297 if ( kk > 0 ) {
298 /*
299 -------------------------------------------
300 load remaining entries in the buffer into B
301 -------------------------------------------
302 */
303 if ( INPMTX_IS_REAL_ENTRIES(A) ) {
304 InpMtx_inputRealTriples(B, kk, ibuf1, ibuf2, dbuf) ;
305 } else if ( INPMTX_IS_COMPLEX_ENTRIES(A) ) {
306 InpMtx_inputComplexTriples(B, kk, ibuf1, ibuf2, dbuf) ;
307 } else {
308 InpMtx_inputTriples(B, kk, ibuf1, ibuf2) ;
309 }
310 }
311 if ( msglvl > 2 ) {
312 fprintf(msgFile, "\n B") ;
313 InpMtx_writeForHumanEye(B, msgFile) ;
314 fflush(msgFile) ;
315 }
316 /*
317 -----------------------------------
318 set B's coordinate type and storage
319 mode to be the same as A's on input
320 -----------------------------------
321 */
322 InpMtx_changeCoordType(B, oldCoordType) ;
323 InpMtx_changeStorageMode(B, oldStorageMode) ;
324 /*
325 -------------------------------------------
326 change back to the old coordinate type of A
327 -------------------------------------------
328 */
329 if ( (oldCoordType = A->coordType) != INPMTX_BY_ROWS ) {
330 InpMtx_changeCoordType(A, oldCoordType) ;
331 }
332 if ( oldStorageMode != A->storageMode ) {
333 InpMtx_changeStorageMode(A, oldStorageMode) ;
334 }
335 /*
336 ------------------------
337 free the working storage
338 ------------------------
339 */
340 IVfree(colmap) ;
341 IVfree(ibuf1) ;
342 IVfree(ibuf2) ;
343 if ( dbuf != NULL ) {
344 DVfree(dbuf) ;
345 }
346 return(1) ; }
347
348 /*--------------------------------------------------------------------*/
349