1 /* copyEntriesToVector.c */
2
3 #include "../A2.h"
4 #include "../../Drand.h"
5
6 /*--------------------------------------------------------------------*/
7 /*
8 ----------------------------------------------------------------
9 purpose -- copy entries to a vector. the portion copied
10 can be a union of the strict lower portion,
11 the diagonal portion, and the strict upper
12 portion. there is one restriction, if the strict
13 lower and strict upper are to be copied, the
14 diagonal will also be copied.
15
16 length -- length of dvec[]
17 dvec[] -- vector to receive matrix entries
18 copyflag -- flag to denote what part of the entries to move
19 A2_STRICT_LOWER --> move strict lower entries
20 A2_LOWER --> move lower entries (includes the diagonal)
21 A2_DIAGONAL --> move diagonal entries
22 A2_UPPER --> move upper entries (includes the diagonal)
23 A2_STRICT_UPPER --> move strict upper entries
24 A2_ALL_ENTRIES --> move all entries
25 storeflag -- flag to denote how to store entries in dvec[]
26 A2_BY_ROWS --> store by rows
27 A2_BY_COLUMNS --> store by columns
28
29 return value -- number of entries copied
30
31 created -- 97jun03, cca, dkw
32 modified -- 98may25, cca
33 ----------------------------------------------------------------
34 */
35 int
A2_copyEntriesToVector(A2 * mtx,int length,double * dvec,int copyflag,int storeflag)36 A2_copyEntriesToVector (
37 A2 *mtx,
38 int length,
39 double *dvec,
40 int copyflag,
41 int storeflag
42 ) {
43 int inc1, inc2, kk, ncol, ndiag, nent, nrow ;
44 /*
45 --------------------------------------------
46 check the input, get dimensions and pointers
47 and check that length is large enough
48 --------------------------------------------
49 */
50 if ( mtx == NULL || length < 0 || dvec == NULL ) {
51 fprintf(stderr,
52 "\n fatal error in A2_copyEntriesToVector(%p,%d,%p,,%d,%d)"
53 "\n bad input\n", mtx, length, dvec, copyflag, storeflag) ;
54 exit(-1) ;
55 }
56 if ( copyflag < 1 || copyflag > 6 ) {
57 fprintf(stderr,
58 "\n fatal error in A2_copyEntriesToVector(%p,%d,%p,%d,%d)"
59 "\n bad input\n"
60 "\n copyflag = %d, must be\n"
61 "\n 1 --> strictly lower entries"
62 "\n 2 --> lower entries"
63 "\n 3 --> diagonal entries"
64 "\n 4 --> strictly upper entries"
65 "\n 5 --> upper entries"
66 "\n 6 --> all entries",
67 mtx, length, dvec, copyflag, storeflag, copyflag) ;
68 exit(-1) ;
69 }
70 if ( storeflag < 0 || storeflag > 1 ) {
71 fprintf(stderr,
72 "\n fatal error in A2_copyEntriesToVector(%p,%d,%p,%d,%d)"
73 "\n bad input\n"
74 "\n storeflag = %d, must be\n"
75 "\n 0 --> store by rows"
76 "\n 1 --> store by columns",
77 mtx, length, dvec, copyflag, storeflag, storeflag) ;
78 exit(-1) ;
79 }
80 nrow = mtx->n1 ;
81 ncol = mtx->n2 ;
82 inc1 = mtx->inc1 ;
83 inc2 = mtx->inc2 ;
84 if ( nrow >= ncol ) {
85 ndiag = ncol ;
86 } else {
87 ndiag = nrow ;
88 }
89 /*
90 ------------------------------------------
91 compute the number of entries to be copied
92 ------------------------------------------
93 */
94 switch ( copyflag ) {
95 case A2_STRICT_LOWER : /* strictly lower entries */
96 nent = (ndiag*(ndiag - 1))/2 + (nrow - ndiag)*ncol ;
97 break ;
98 case A2_LOWER : /* lower entries */
99 nent = (ndiag*(ndiag + 1))/2 + (nrow - ndiag)*ncol ;
100 break ;
101 case A2_DIAGONAL : /* diagonal entries */
102 nent = ndiag ;
103 break ;
104 case A2_UPPER : /* upper entries */
105 nent = (ndiag*(ndiag + 1))/2 + (ncol - ndiag)*nrow ;
106 break ;
107 case A2_STRICT_UPPER : /* strictly upper entries */
108 nent = (ndiag*(ndiag - 1))/2 + (ncol - ndiag)*nrow ;
109 break ;
110 case A2_ALL_ENTRIES : /* all entries */
111 nent = nrow*ncol ;
112 break ;
113 default :
114 break ;
115 }
116 if ( nent > length ) {
117 fprintf(stderr,
118 "\n fatal error in A2_copyEntriesToVector(%p,%d,%p,%d,%d)"
119 "\n nent = %d, buffer length = %d",
120 mtx, length, dvec, copyflag, storeflag, nent, length) ;
121 exit(-1) ;
122 }
123 /*
124 --------------------------------------------
125 make life simple, unit stride through dvec[]
126 --------------------------------------------
127 */
128 kk = 0 ;
129 if ( storeflag == A2_BY_ROWS ) {
130 int irow, jstart, jcol, jend, jj ;
131 double *row ;
132 /*
133 --------------
134 loop over rows
135 --------------
136 */
137 switch ( copyflag ) {
138 case A2_STRICT_LOWER :
139 /*
140 ----------------------
141 strictly lower entries
142 ----------------------
143 */
144 if ( A2_IS_REAL(mtx) ) {
145 for ( irow = 0, row = mtx->entries, kk = 0 ;
146 irow < nrow ;
147 irow++, row += inc1 ) {
148 jstart = 0 ;
149 jend = (irow < ndiag) ? irow - 1 : ndiag - 1 ;
150 for ( jcol = jstart, jj = jcol*inc2 ;
151 jcol <= jend ;
152 jcol++, jj += inc2, kk++ ) {
153 dvec[kk] = row[jj] ;
154 }
155 }
156 } else if ( A2_IS_COMPLEX(mtx) ) {
157 for ( irow = 0, row = mtx->entries, kk = 0 ;
158 irow < nrow ;
159 irow++, row += 2*inc1 ) {
160 jstart = 0 ;
161 jend = (irow < ndiag) ? irow - 1 : ndiag - 1 ;
162 for ( jcol = jstart, jj = jcol*inc2 ;
163 jcol <= jend ;
164 jcol++, jj += inc2, kk++ ) {
165 dvec[2*kk] = row[2*jj] ;
166 dvec[2*kk+1] = row[2*jj+1] ;
167 }
168 }
169 }
170 break ;
171 case A2_LOWER :
172 /*
173 ------------------------------------
174 lower entries including the diagonal
175 ------------------------------------
176 */
177 if ( A2_IS_REAL(mtx) ) {
178 for ( irow = 0, row = mtx->entries, kk = 0 ;
179 irow < nrow ;
180 irow++, row += inc1 ) {
181 jstart = 0 ;
182 jend = (irow < ndiag) ? irow : ndiag - 1 ;
183 for ( jcol = jstart, jj = jcol*inc2 ;
184 jcol <= jend ;
185 jcol++, jj += inc2, kk++ ) {
186 dvec[kk] = row[jj] ;
187 }
188 }
189 } else if ( A2_IS_COMPLEX(mtx) ) {
190 for ( irow = 0, row = mtx->entries, kk = 0 ;
191 irow < nrow ;
192 irow++, row += 2*inc1 ) {
193 jstart = 0 ;
194 jend = (irow < ndiag) ? irow : ndiag - 1 ;
195 for ( jcol = jstart, jj = jcol*inc2 ;
196 jcol <= jend ;
197 jcol++, jj += inc2, kk++ ) {
198 dvec[2*kk] = row[2*jj] ;
199 dvec[2*kk+1] = row[2*jj+1] ;
200 }
201 }
202 }
203 break ;
204 case A2_DIAGONAL :
205 /*
206 -----------------
207 just the diagonal
208 -----------------
209 */
210 if ( A2_IS_REAL(mtx) ) {
211 for ( irow = 0, row = mtx->entries, kk = 0 ;
212 irow < ndiag ;
213 irow++, row += inc1, kk++ ) {
214 dvec[kk] = row[irow*inc2] ;
215 }
216 } else if ( A2_IS_COMPLEX(mtx) ) {
217 for ( irow = 0, row = mtx->entries, kk = 0 ;
218 irow < ndiag ;
219 irow++, row += 2*inc1, kk++ ) {
220 dvec[2*kk] = row[2*irow*inc2] ;
221 dvec[2*kk+1] = row[2*irow*inc2+1] ;
222 }
223 }
224 break ;
225 case A2_UPPER :
226 /*
227 ------------------------------------
228 upper entries including the diagonal
229 ------------------------------------
230 */
231 if ( A2_IS_REAL(mtx) ) {
232 for ( irow = 0, row = mtx->entries, kk = 0 ;
233 irow < nrow ;
234 irow++, row += inc1 ) {
235 jstart = irow ;
236 jend = ncol - 1 ;
237 for ( jcol = jstart, jj = jcol*inc2 ;
238 jcol <= jend ;
239 jcol++, jj += inc2, kk++ ) {
240 dvec[kk] = row[jj] ;
241 }
242 }
243 } else if ( A2_IS_COMPLEX(mtx) ) {
244 for ( irow = 0, row = mtx->entries, kk = 0 ;
245 irow < nrow ;
246 irow++, row += 2*inc1 ) {
247 jstart = irow ;
248 jend = ncol - 1 ;
249 for ( jcol = jstart, jj = jcol*inc2 ;
250 jcol <= jend ;
251 jcol++, jj += inc2, kk++ ) {
252 dvec[2*kk] = row[2*jj] ;
253 dvec[2*kk+1] = row[2*jj+1] ;
254 }
255 }
256 }
257 break ;
258 case A2_STRICT_UPPER :
259 /*
260 --------------------------
261 strictly the upper entries
262 --------------------------
263 */
264 if ( A2_IS_REAL(mtx) ) {
265 for ( irow = 0, row = mtx->entries, kk = 0 ;
266 irow < nrow ;
267 irow++, row += inc1 ) {
268 jstart = irow + 1 ;
269 jend = ncol - 1 ;
270 for ( jcol = jstart, jj = jcol*inc2 ;
271 jcol <= jend ;
272 jcol++, jj += inc2, kk++ ) {
273 dvec[kk] = row[jj] ;
274 }
275 }
276 } else if ( A2_IS_COMPLEX(mtx) ) {
277 for ( irow = 0, row = mtx->entries, kk = 0 ;
278 irow < nrow ;
279 irow++, row += 2*inc1 ) {
280 jstart = irow + 1 ;
281 jend = ncol - 1 ;
282 for ( jcol = jstart, jj = jcol*inc2 ;
283 jcol <= jend ;
284 jcol++, jj += inc2, kk++ ) {
285 dvec[2*kk] = row[2*jj] ;
286 dvec[2*kk+1] = row[2*jj+1] ;
287 }
288 }
289 }
290 break ;
291 case A2_ALL_ENTRIES :
292 /*
293 -----------
294 all entries
295 -----------
296 */
297 if ( A2_IS_REAL(mtx) ) {
298 for ( irow = 0, row = mtx->entries, kk = 0 ;
299 irow < nrow ;
300 irow++, row += inc1 ) {
301 jstart = 0 ;
302 jend = ncol - 1 ;
303 for ( jcol = jstart, jj = 0 ;
304 jcol <= jend ;
305 jcol++, jj += inc2, kk++ ) {
306 dvec[kk] = row[jj] ;
307 }
308 }
309 } else if ( A2_IS_COMPLEX(mtx) ) {
310 for ( irow = 0, row = mtx->entries, kk = 0 ;
311 irow < nrow ;
312 irow++, row += 2*inc1 ) {
313 jstart = 0 ;
314 jend = ncol - 1 ;
315 for ( jcol = jstart, jj = 0 ;
316 jcol <= jend ;
317 jcol++, jj += inc2, kk++ ) {
318 dvec[2*kk] = row[2*jj] ;
319 dvec[2*kk+1] = row[2*jj+1] ;
320 }
321 }
322 }
323 break ;
324 default :
325 break ;
326 }
327 } else {
328 int iend, ii, irow, istart, jcol ;
329 double *col ;
330 /*
331 -----------------
332 loop over columns
333 -----------------
334 */
335 kk = 0 ;
336 switch ( copyflag ) {
337 case A2_STRICT_LOWER :
338 /*
339 ----------------------
340 strictly lower entries
341 ----------------------
342 */
343 if ( A2_IS_REAL(mtx) ) {
344 for ( jcol = 0, col = mtx->entries, kk = 0 ;
345 jcol < ncol ;
346 jcol++, col += inc2 ) {
347 istart = jcol + 1 ;
348 iend = nrow - 1 ;
349 for ( irow = istart, ii = irow*inc1 ;
350 irow <= iend ;
351 irow++, ii += inc1, kk++ ) {
352 dvec[kk] = col[ii] ;
353 }
354 }
355 } else if ( A2_IS_COMPLEX(mtx) ) {
356 for ( jcol = 0, col = mtx->entries, kk = 0 ;
357 jcol < ncol ;
358 jcol++, col += 2*inc2 ) {
359 istart = jcol + 1 ;
360 iend = nrow - 1 ;
361 for ( irow = istart, ii = irow*inc1 ;
362 irow <= iend ;
363 irow++, ii += inc1, kk++ ) {
364 dvec[2*kk] = col[2*ii] ;
365 dvec[2*kk+1] = col[2*ii+1] ;
366 }
367 }
368 }
369 break ;
370 case A2_LOWER :
371 /*
372 ------------------------------------
373 lower entries including the diagonal
374 ------------------------------------
375 */
376 if ( A2_IS_REAL(mtx) ) {
377 for ( jcol = 0, col = mtx->entries, kk = 0 ;
378 jcol < ncol ;
379 jcol++, col += inc2 ) {
380 istart = jcol ;
381 iend = nrow - 1 ;
382 for ( irow = istart, ii = irow*inc1 ;
383 irow <= iend ;
384 irow++, ii += inc1, kk++ ) {
385 dvec[kk] = col[ii] ;
386 }
387 }
388 } else if ( A2_IS_COMPLEX(mtx) ) {
389 for ( jcol = 0, col = mtx->entries, kk = 0 ;
390 jcol < ncol ;
391 jcol++, col += 2*inc2 ) {
392 istart = jcol ;
393 iend = nrow - 1 ;
394 for ( irow = istart, ii = irow*inc1 ;
395 irow <= iend ;
396 irow++, ii += inc1, kk++ ) {
397 dvec[2*kk] = col[2*ii] ;
398 dvec[2*kk+1] = col[2*ii+1] ;
399 }
400 }
401 }
402 break ;
403 case A2_DIAGONAL :
404 /*
405 -----------------
406 just the diagonal
407 -----------------
408 */
409 if ( A2_IS_REAL(mtx) ) {
410 for ( jcol = 0, col = mtx->entries, kk = 0 ;
411 jcol < ndiag ;
412 jcol++, col += inc2, kk++ ) {
413 dvec[kk] = col[jcol*inc1] ;
414 }
415 } else if ( A2_IS_COMPLEX(mtx) ) {
416 for ( jcol = 0, col = mtx->entries, kk = 0 ;
417 jcol < ndiag ;
418 jcol++, col += 2*inc2, kk++ ) {
419 dvec[2*kk] = col[2*jcol*inc1] ;
420 dvec[2*kk+1] = col[2*jcol*inc1+1] ;
421 }
422 }
423 break ;
424 case A2_UPPER :
425 /*
426 ------------------------------------
427 upper entries including the diagonal
428 ------------------------------------
429 */
430 if ( A2_IS_REAL(mtx) ) {
431 for ( jcol = 0, col = mtx->entries, kk = 0 ;
432 jcol < ncol ;
433 jcol++, col += inc2 ) {
434 istart = 0 ;
435 iend = (jcol < ndiag) ? jcol : ndiag - 1 ;
436 for ( irow = istart, ii = irow*inc1 ;
437 irow <= iend ;
438 irow++, ii += inc1, kk++ ) {
439 dvec[kk] = col[ii] ;
440 }
441 }
442 } else if ( A2_IS_COMPLEX(mtx) ) {
443 for ( jcol = 0, col = mtx->entries, kk = 0 ;
444 jcol < ncol ;
445 jcol++, col += 2*inc2 ) {
446 istart = 0 ;
447 iend = (jcol < ndiag) ? jcol : ndiag - 1 ;
448 for ( irow = istart, ii = irow*inc1 ;
449 irow <= iend ;
450 irow++, ii += inc1, kk++ ) {
451 dvec[2*kk] = col[2*ii] ;
452 dvec[2*kk+1] = col[2*ii+1] ;
453 }
454 }
455 }
456 break ;
457 case A2_STRICT_UPPER :
458 /*
459 --------------------------
460 strictly the upper entries
461 --------------------------
462 */
463 if ( A2_IS_REAL(mtx) ) {
464 for ( jcol = 0, col = mtx->entries, kk = 0 ;
465 jcol < ncol ;
466 jcol++, col += inc2 ) {
467 istart = 0 ;
468 iend = (jcol < ndiag) ? jcol - 1 : ndiag - 1 ;
469 for ( irow = istart, ii = irow*inc1 ;
470 irow <= iend ;
471 irow++, ii += inc1, kk++ ) {
472 dvec[kk] = col[ii] ;
473 }
474 }
475 } else if ( A2_IS_COMPLEX(mtx) ) {
476 for ( jcol = 0, col = mtx->entries, kk = 0 ;
477 jcol < ncol ;
478 jcol++, col += 2*inc2 ) {
479 istart = 0 ;
480 iend = (jcol < ndiag) ? jcol - 1 : ndiag - 1 ;
481 for ( irow = istart, ii = irow*inc1 ;
482 irow <= iend ;
483 irow++, ii += inc1, kk++ ) {
484 dvec[2*kk] = col[2*ii] ;
485 dvec[2*kk+1] = col[2*ii+1] ;
486 }
487 }
488 }
489 break ;
490 case A2_ALL_ENTRIES :
491 /*
492 -----------
493 all entries
494 -----------
495 */
496 if ( A2_IS_REAL(mtx) ) {
497 for ( jcol = 0, col = mtx->entries, kk = 0 ;
498 jcol < ncol ;
499 jcol++, col += inc2 ) {
500 istart = 0 ;
501 iend = nrow - 1 ;
502 for ( irow = istart, ii = 0 ;
503 irow <= iend ;
504 irow++, ii += inc1, kk++ ) {
505 dvec[kk] = col[ii] ;
506 }
507 }
508 } else if ( A2_IS_COMPLEX(mtx) ) {
509 for ( jcol = 0, col = mtx->entries, kk = 0 ;
510 jcol < ncol ;
511 jcol++, col += 2*inc2 ) {
512 istart = 0 ;
513 iend = nrow - 1 ;
514 for ( irow = istart, ii = 0 ;
515 irow <= iend ;
516 irow++, ii += inc1, kk++ ) {
517 dvec[2*kk] = col[2*ii] ;
518 dvec[2*kk+1] = col[2*ii+1] ;
519 }
520 }
521 }
522 break ;
523 default :
524 break ;
525 }
526 }
527 return(kk) ; }
528
529 /*--------------------------------------------------------------------*/
530