1 /* util.c */
2
3 #include "../A2.h"
4 #include "../../Drand.h"
5
6 /*--------------------------------------------------------------------*/
7 /*
8 ----------------------------------------------
9 return the number of bytes taken by the object
10
11 created -- 98may01, cca
12 ----------------------------------------------
13 */
14 int
A2_sizeOf(A2 * mtx)15 A2_sizeOf (
16 A2 *mtx
17 ) {
18 int nbytes ;
19 /*
20 ---------------
21 check the input
22 ---------------
23 */
24 if ( mtx == NULL ) {
25 fprintf(stderr, "\n fatal error in A2_sizeOf(%p)"
26 "\n bad input\n", mtx) ;
27 exit(-1) ;
28 }
29 if ( ! (A2_IS_REAL(mtx) || A2_IS_COMPLEX(mtx)) ) {
30 fprintf(stderr, "\n fatal error in A2_sizeOf(%p)"
31 "\n bad type %d, must be SPOOLES_REAL or SPOOLES_COMPLEX\n",
32 mtx, mtx->type) ;
33 exit(-1) ;
34 }
35 if ( A2_IS_REAL(mtx) ) {
36 nbytes = sizeof(struct _A2) + mtx->nowned*sizeof(double) ;
37 } else if ( A2_IS_COMPLEX(mtx) ) {
38 nbytes = sizeof(struct _A2) + 2*mtx->nowned*sizeof(double) ;
39 }
40 return(nbytes) ; }
41
42 /*--------------------------------------------------------------------*/
43 /*
44 ---------------------------------------------------------------
45 shift the base of the entries and adjust dimensions
46
47 mtx(0:n1-rowoff-1,0:n2-coloff-1) = mtx(rowoff:n1-1,coloff:n2-1)
48
49 created -- 98may01, cca
50 ---------------------------------------------------------------
51 */
52 void
A2_shiftBase(A2 * mtx,int rowoff,int coloff)53 A2_shiftBase (
54 A2 *mtx,
55 int rowoff,
56 int coloff
57 ) {
58 /*
59 ---------------
60 check the input
61 ---------------
62 */
63 if ( mtx == NULL ) {
64 fprintf(stderr, "\n fatal error in A2_shiftbase(%p,%d,%d)"
65 "\n bad input\n", mtx, rowoff, coloff) ;
66 exit(-1) ;
67 }
68 if ( ! (A2_IS_REAL(mtx) || A2_IS_COMPLEX(mtx)) ) {
69 fprintf(stderr, "\n fatal error in A2_shiftBase(%p,%d,%d)"
70 "\n bad type %d, must be SPOOLES_REAL or SPOOLES_COMPLEX\n",
71 mtx, rowoff, coloff, mtx->type) ;
72 exit(-1) ;
73 }
74 mtx->n1 -= rowoff ;
75 mtx->n2 -= coloff ;
76 if ( A2_IS_REAL(mtx) ) {
77 mtx->entries += rowoff*mtx->inc1 + coloff*mtx->inc2 ;
78 } else if ( A2_IS_COMPLEX(mtx) ) {
79 mtx->entries += 2*(rowoff*mtx->inc1 + coloff*mtx->inc2) ;
80 }
81 return ; }
82
83 /*--------------------------------------------------------------------*/
84 /*
85 --------------------------------------------------------------
86 returns 1 if the storage is row major, otherwise returns zero.
87
88 created -- 98may01, cca
89 --------------------------------------------------------------
90 */
91 int
A2_rowMajor(A2 * mtx)92 A2_rowMajor (
93 A2 *mtx
94 ) {
95 /*
96 ---------------
97 check the input
98 ---------------
99 */
100 if ( mtx == NULL ) {
101 fprintf(stderr, "\n fatal error in A2_rowMajor(%p)"
102 "\n bad input\n", mtx) ;
103 exit(-1) ;
104 }
105 if ( ! (A2_IS_REAL(mtx) || A2_IS_COMPLEX(mtx)) ) {
106 fprintf(stderr, "\n fatal error in A2_rowMajor(%p)"
107 "\n bad type %d, must be SPOOLES_REAL or SPOOLES_COMPLEX\n",
108 mtx, mtx->type) ;
109 exit(-1) ;
110 }
111 if ( mtx->inc2 == 1 ) {
112 return(1) ;
113 } else {
114 return(0) ;
115 } }
116
117 /*--------------------------------------------------------------------*/
118 /*
119 -----------------------------------------------------------------
120 returns 1 if the storage is column major, otherwise returns zero.
121
122 created -- 98may01, cca
123 -----------------------------------------------------------------
124 */
125 int
A2_columnMajor(A2 * mtx)126 A2_columnMajor (
127 A2 *mtx
128 ) {
129 /*
130 ---------------
131 check the input
132 ---------------
133 */
134 if ( mtx == NULL ) {
135 fprintf(stderr, "\n fatal error in A2_columnMajor(%p)"
136 "\n bad input\n", mtx) ;
137 exit(-1) ;
138 }
139 if ( ! (A2_IS_REAL(mtx) || A2_IS_COMPLEX(mtx)) ) {
140 fprintf(stderr, "\n fatal error in A2_columnMajor(%p)"
141 "\n bad type %d, must be SPOOLES_REAL or SPOOLES_COMPLEX\n",
142 mtx, mtx->type) ;
143 exit(-1) ;
144 }
145 if ( mtx->inc1 == 1 ) {
146 return(1) ;
147 } else {
148 return(0) ;
149 } }
150
151 /*--------------------------------------------------------------------*/
152 /*
153 -----------------------
154 transpose the matrix
155
156 created -- 98may01, cca
157 -----------------------
158 */
159 void
A2_transpose(A2 * mtx)160 A2_transpose (
161 A2 *mtx
162 ) {
163 int inc1, n1 ;
164 /*
165 ---------------
166 check the input
167 ---------------
168 */
169 if ( mtx == NULL ) {
170 fprintf(stderr, "\n fatal error in A2_transpose(%p)"
171 "\n bad input\n", mtx) ;
172 exit(-1) ;
173 }
174 if ( ! (A2_IS_REAL(mtx) || A2_IS_COMPLEX(mtx)) ) {
175 fprintf(stderr, "\n fatal error in A2_transpose(%p)"
176 "\n bad type %d, must be SPOOLES_REAL or SPOOLES_COMPLEX\n",
177 mtx, mtx->type) ;
178 exit(-1) ;
179 }
180 n1 = mtx->n1 ;
181 mtx->n1 = mtx->n2 ;
182 mtx->n2 = n1 ;
183 inc1 = mtx->inc1 ;
184 mtx->inc1 = mtx->inc2 ;
185 mtx->inc2 = inc1 ;
186
187 return ; }
188
189 /*--------------------------------------------------------------------*/
190 /*
191 ----------------------------
192 extract row[*] = mtx(irow,*)
193
194 created -- 98may01, cca
195 ----------------------------
196 */
197 void
A2_extractRow(A2 * mtx,double row[],int irow)198 A2_extractRow (
199 A2 *mtx,
200 double row[],
201 int irow
202 ) {
203 double *entries ;
204 int inc2, j, k, n2 ;
205 /*
206 ---------------
207 check the input
208 ---------------
209 */
210 if ( mtx == NULL || row == NULL || mtx->entries == NULL
211 || irow < 0 || irow >= mtx->n1 ) {
212 fprintf(stderr, "\n fatal error in A2_extractRow(%p,%p,%d)"
213 "\n bad input\n", mtx, row, irow) ;
214 exit(-1) ;
215 }
216 if ( ! (A2_IS_REAL(mtx) || A2_IS_COMPLEX(mtx)) ) {
217 fprintf(stderr, "\n fatal error in A2_extractRow(%p,%p,%d)"
218 "\n bad type %d, must be SPOOLES_REAL or SPOOLES_COMPLEX\n",
219 mtx, row, irow, mtx->type) ;
220 exit(-1) ;
221 }
222 k = irow * mtx->inc1 ;
223 n2 = mtx->n2 ;
224 inc2 = mtx->inc2 ;
225 entries = mtx->entries ;
226 if ( A2_IS_REAL(mtx) ) {
227 for ( j = 0 ; j < n2 ; j++, k += inc2 ) {
228 row[j] = entries[k] ;
229 }
230 } else if ( A2_IS_COMPLEX(mtx) ) {
231 for ( j = 0 ; j < n2 ; j++, k += inc2 ) {
232 row[2*j] = entries[2*k] ;
233 row[2*j+1] = entries[2*k+1] ;
234 }
235 }
236 return ; }
237
238 /*--------------------------------------------------------------------*/
239 /*
240 ----------------------------
241 extract col[*] = mtx(*,jcol)
242
243 created -- 98may01, cca
244 ----------------------------
245 */
246 void
A2_extractColumn(A2 * mtx,double col[],int jcol)247 A2_extractColumn (
248 A2 *mtx,
249 double col[],
250 int jcol
251 ) {
252 double *entries ;
253 int i, inc1, k, n1 ;
254 /*
255 ---------------
256 check the input
257 ---------------
258 */
259 if ( mtx == NULL || col == NULL || mtx->entries == NULL
260 || jcol < 0 || jcol >= mtx->n2 ) {
261 fprintf(stderr, "\n fatal error in A2_extractColumn(%p,%p,%d)"
262 "\n bad input\n", mtx, col, jcol) ;
263 exit(-1) ;
264 }
265 if ( ! (A2_IS_REAL(mtx) || A2_IS_COMPLEX(mtx)) ) {
266 fprintf(stderr, "\n fatal error in A2_extractColumn(%p,%p,%d)"
267 "\n bad type %d, must be SPOOLES_REAL or SPOOLES_COMPLEX\n",
268 mtx, col, jcol, mtx->type) ;
269 exit(-1) ;
270 }
271 k = jcol * mtx->inc2 ;
272 n1 = mtx->n1 ;
273 inc1 = mtx->inc1 ;
274 entries = mtx->entries ;
275 if ( A2_IS_REAL(mtx) ) {
276 for ( i = 0 ; i < n1 ; i++, k += inc1 ) {
277 col[i] = entries[k] ;
278 }
279 } else if ( A2_IS_COMPLEX(mtx) ) {
280 for ( i = 0 ; i < n1 ; i++, k += inc1 ) {
281 col[2*i] = entries[2*k] ;
282 col[2*i+1] = entries[2*k+1] ;
283 }
284 }
285 return ; }
286
287 /*--------------------------------------------------------------------*/
288 /*
289 -----------------------
290 set mtx(irow,*) = y[*]
291
292 created -- 98may01, cca
293 -----------------------
294 */
295 void
A2_setRow(A2 * mtx,double row[],int irow)296 A2_setRow (
297 A2 *mtx,
298 double row[],
299 int irow
300 ) {
301 double *entries ;
302 int inc2, j, k, n2 ;
303 /*
304 ---------------
305 check the input
306 ---------------
307 */
308 if ( mtx == NULL || row == NULL || irow < 0 || irow >= mtx->n1 ) {
309 fprintf(stderr, "\n fatal error in A2_setRow(%p,%p,%d)"
310 "\n bad input\n", mtx, row, irow) ;
311 exit(-1) ;
312 }
313 if ( ! (A2_IS_REAL(mtx) || A2_IS_COMPLEX(mtx)) ) {
314 fprintf(stderr, "\n fatal error in A2_setRow(%p,%p,%d)"
315 "\n bad type %d, must be SPOOLES_REAL or SPOOLES_COMPLEX\n",
316 mtx, row, irow, mtx->type) ;
317 exit(-1) ;
318 }
319 k = irow * mtx->inc1 ;
320 n2 = mtx->n2 ;
321 inc2 = mtx->inc2 ;
322 entries = mtx->entries ;
323 if ( A2_IS_REAL(mtx) ) {
324 for ( j = 0 ; j < n2 ; j++, k += inc2 ) {
325 entries[k] = row[j] ;
326 }
327 } else if ( A2_IS_COMPLEX(mtx) ) {
328 for ( j = 0 ; j < n2 ; j++, k += inc2 ) {
329 entries[2*k] = row[2*j] ;
330 entries[2*k+1] = row[2*j+1] ;
331 }
332 }
333 return ; }
334
335 /*--------------------------------------------------------------------*/
336 /*
337 -----------------------
338 set mtx(*,jcol) = y[*]
339
340 created -- 98may01, cca
341 -----------------------
342 */
343 void
A2_setColumn(A2 * mtx,double col[],int jcol)344 A2_setColumn (
345 A2 *mtx,
346 double col[],
347 int jcol
348 ) {
349 double *entries ;
350 int inc1, i, k, n1 ;
351 /*
352 ---------------
353 check the input
354 ---------------
355 */
356 if ( mtx == NULL || col == NULL || jcol < 0 || jcol >= mtx->n2 ) {
357 fprintf(stderr, "\n fatal error in A2_setColumn(%p,%p,%d)"
358 "\n bad input\n", mtx, col, jcol) ;
359 exit(-1) ;
360 }
361 if ( ! (A2_IS_REAL(mtx) || A2_IS_COMPLEX(mtx)) ) {
362 fprintf(stderr, "\n fatal error in A2_setColumn(%p,%p,%d)"
363 "\n bad type %d, must be SPOOLES_REAL or SPOOLES_COMPLEX\n",
364 mtx, col, jcol, mtx->type) ;
365 exit(-1) ;
366 }
367 k = jcol * mtx->inc2 ;
368 n1 = mtx->n1 ;
369 inc1 = mtx->inc1 ;
370 entries = mtx->entries ;
371 if ( A2_IS_REAL(mtx) ) {
372 for ( i = 0 ; i < n1 ; i++, k += inc1 ) {
373 entries[k] = col[i] ;
374 }
375 } else if ( A2_IS_COMPLEX(mtx) ) {
376 for ( i = 0 ; i < n1 ; i++, k += inc1 ) {
377 entries[2*k] = col[2*i] ;
378 entries[2*k+1] = col[2*i+1] ;
379 }
380 }
381 return ; }
382
383 /*--------------------------------------------------------------------*/
384 /*
385 ----------------------------
386 extract row[*] = mtx(irow,*)
387
388 created -- 98may01, cca
389 ----------------------------
390 */
391 void
A2_extractRowDV(A2 * mtx,DV * rowDV,int irow)392 A2_extractRowDV (
393 A2 *mtx,
394 DV *rowDV,
395 int irow
396 ) {
397 double *entries, *row ;
398 int inc2, j, k, n2 ;
399 /*
400 ---------------
401 check the input
402 ---------------
403 */
404 if ( mtx == NULL || rowDV == NULL || mtx->entries == NULL
405 || irow < 0 || irow >= mtx->n1 ) {
406 fprintf(stderr, "\n fatal error in A2_extractRowDV(%p,%p,%d)"
407 "\n bad input\n", mtx, rowDV, irow) ;
408 exit(-1) ;
409 }
410 if ( ! A2_IS_REAL(mtx) ) {
411 fprintf(stderr, "\n fatal error in A2_extractRowDV(%p,%p,%d)"
412 "\n bad type %d, must be SPOOLES_REAL\n",
413 mtx, rowDV, irow, mtx->type) ;
414 exit(-1) ;
415 }
416 if ( DV_size(rowDV) < (n2 = mtx->n2) ) {
417 DV_setSize(rowDV, n2) ;
418 }
419 row = DV_entries(rowDV) ;
420 k = irow * mtx->inc1 ;
421 inc2 = mtx->inc2 ;
422 entries = mtx->entries ;
423 for ( j = 0 ; j < n2 ; j++, k += inc2 ) {
424 row[j] = entries[k] ;
425 }
426 return ; }
427
428 /*--------------------------------------------------------------------*/
429 /*
430 ----------------------------
431 extract row[*] = mtx(irow,*)
432
433 created -- 98may01, cca
434 ----------------------------
435 */
436 void
A2_extractRowZV(A2 * mtx,ZV * rowZV,int irow)437 A2_extractRowZV (
438 A2 *mtx,
439 ZV *rowZV,
440 int irow
441 ) {
442 double *entries, *row ;
443 int inc2, j, k, n2 ;
444 /*
445 ---------------
446 check the input
447 ---------------
448 */
449 if ( mtx == NULL || rowZV == NULL || mtx->entries == NULL
450 || irow < 0 || irow >= mtx->n1 ) {
451 fprintf(stderr, "\n fatal error in A2_extractRowZV(%p,%p,%d)"
452 "\n bad input\n", mtx, rowZV, irow) ;
453 exit(-1) ;
454 }
455 if ( ! A2_IS_COMPLEX(mtx) ) {
456 fprintf(stderr, "\n fatal error in A2_extractRowZV(%p,%p,%d)"
457 "\n bad type %d, must be SPOOLES_COMPLEX\n",
458 mtx, rowZV, irow, mtx->type) ;
459 exit(-1) ;
460 }
461 if ( ZV_size(rowZV) < (n2 = mtx->n2) ) {
462 ZV_setSize(rowZV, n2) ;
463 }
464 row = ZV_entries(rowZV) ;
465 k = irow * mtx->inc1 ;
466 inc2 = mtx->inc2 ;
467 entries = mtx->entries ;
468 for ( j = 0 ; j < n2 ; j++, k += inc2 ) {
469 row[2*j] = entries[2*k] ;
470 row[2*j+1] = entries[2*k+1] ;
471 }
472 return ; }
473
474 /*--------------------------------------------------------------------*/
475 /*
476 ----------------------------
477 extract col[*] = mtx(*,jcol)
478
479 created -- 98may01, cca
480 ----------------------------
481 */
482 void
A2_extractColumnDV(A2 * mtx,DV * colDV,int jcol)483 A2_extractColumnDV (
484 A2 *mtx,
485 DV *colDV,
486 int jcol
487 ) {
488 double *entries, *col ;
489 int i, inc1, k, n1 ;
490 /*
491 ---------------
492 check the input
493 ---------------
494 */
495 if ( mtx == NULL || colDV == NULL || mtx->entries == NULL
496 || jcol < 0 || jcol >= mtx->n2 ) {
497 fprintf(stderr, "\n fatal error in A2_extractColumnDV(%p,%p,%d)"
498 "\n bad input\n", mtx, colDV, jcol) ;
499 exit(-1) ;
500 }
501 if ( ! A2_IS_REAL(mtx) ) {
502 fprintf(stderr, "\n fatal error in A2_extractColumnDV(%p,%p,%d)"
503 "\n bad type %d, must be SPOOLES_REAL\n",
504 mtx, colDV, jcol, mtx->type) ;
505 exit(-1) ;
506 }
507 if ( DV_size(colDV) < (n1 = mtx->n1) ) {
508 DV_setSize(colDV, n1) ;
509 }
510 col = DV_entries(colDV) ;
511 k = jcol * mtx->inc2 ;
512 inc1 = mtx->inc1 ;
513 entries = mtx->entries ;
514 for ( i = 0 ; i < n1 ; i++, k += inc1 ) {
515 col[i] = entries[k] ;
516 }
517 return ; }
518
519 /*--------------------------------------------------------------------*/
520 /*
521 ----------------------------
522 extract col[*] = mtx(*,jcol)
523
524 created -- 98may01, cca
525 ----------------------------
526 */
527 void
A2_extractColumnZV(A2 * mtx,ZV * colZV,int jcol)528 A2_extractColumnZV (
529 A2 *mtx,
530 ZV *colZV,
531 int jcol
532 ) {
533 double *entries, *col ;
534 int i, inc1, k, n1 ;
535 /*
536 ---------------
537 check the input
538 ---------------
539 */
540 if ( mtx == NULL || colZV == NULL || mtx->entries == NULL
541 || jcol < 0 || jcol >= mtx->n2 ) {
542 fprintf(stderr, "\n fatal error in A2_extractColumnZV(%p,%p,%d)"
543 "\n bad input\n", mtx, colZV, jcol) ;
544 exit(-1) ;
545 }
546 if ( ! A2_IS_COMPLEX(mtx) ) {
547 fprintf(stderr, "\n fatal error in A2_extractColumnZV(%p,%p,%d)"
548 "\n bad type %d, must be SPOOLES_COMPLEX\n",
549 mtx, colZV, jcol, mtx->type) ;
550 exit(-1) ;
551 }
552 if ( ZV_size(colZV) < (n1 = mtx->n1) ) {
553 ZV_setSize(colZV, n1) ;
554 }
555 col = ZV_entries(colZV) ;
556 k = jcol * mtx->inc2 ;
557 inc1 = mtx->inc1 ;
558 entries = mtx->entries ;
559 for ( i = 0 ; i < n1 ; i++, k += inc1 ) {
560 col[2*i] = entries[2*k] ;
561 col[2*i+1] = entries[2*k+1] ;
562 }
563 return ; }
564
565 /*--------------------------------------------------------------------*/
566 /*
567 -----------------------
568 set mtx(irow,*) = y[*]
569
570 created -- 98may01, cca
571 -----------------------
572 */
573 void
A2_setRowDV(A2 * mtx,DV * rowDV,int irow)574 A2_setRowDV (
575 A2 *mtx,
576 DV *rowDV,
577 int irow
578 ) {
579 double *entries, *row ;
580 int inc2, j, k, n2 ;
581 /*
582 ---------------
583 check the input
584 ---------------
585 */
586 if ( mtx == NULL || rowDV == NULL || DV_size(rowDV) != (n2 = mtx->n2)
587 || irow < 0 || irow >= mtx->n1 ) {
588 fprintf(stderr, "\n fatal error in A2_setRowDV(%p,%p,%d)"
589 "\n bad input\n", mtx, rowDV, irow) ;
590 exit(-1) ;
591 }
592 if ( ! A2_IS_REAL(mtx) ) {
593 fprintf(stderr, "\n fatal error in A2_setRowDV(%p,%p,%d)"
594 "\n bad type %d, must be SPOOLES_REAL\n",
595 mtx, rowDV, irow, mtx->type) ;
596 exit(-1) ;
597 }
598 k = irow * mtx->inc1 ;
599 inc2 = mtx->inc2 ;
600 entries = mtx->entries ;
601 row = DV_entries(rowDV) ;
602 for ( j = 0 ; j < n2 ; j++, k += inc2 ) {
603 entries[k] = row[j] ;
604 }
605 return ; }
606
607 /*--------------------------------------------------------------------*/
608 /*
609 -----------------------
610 set mtx(irow,*) = y[*]
611
612 created -- 98may01, cca
613 -----------------------
614 */
615 void
A2_setRowZV(A2 * mtx,ZV * rowZV,int irow)616 A2_setRowZV (
617 A2 *mtx,
618 ZV *rowZV,
619 int irow
620 ) {
621 double *entries, *row ;
622 int inc2, j, k, n2 ;
623 /*
624 ---------------
625 check the input
626 ---------------
627 */
628 if ( mtx == NULL || rowZV == NULL || ZV_size(rowZV) != (n2 = mtx->n2)
629 || irow < 0 || irow >= mtx->n1 ) {
630 fprintf(stderr, "\n fatal error in A2_setRowZV(%p,%p,%d)"
631 "\n bad input\n", mtx, rowZV, irow) ;
632 exit(-1) ;
633 }
634 if ( ! A2_IS_COMPLEX(mtx) ) {
635 fprintf(stderr, "\n fatal error in A2_setRowZV(%p,%p,%d)"
636 "\n bad type %d, must be SPOOLES_COMPLEX\n",
637 mtx, rowZV, irow, mtx->type) ;
638 exit(-1) ;
639 }
640 k = irow * mtx->inc1 ;
641 inc2 = mtx->inc2 ;
642 entries = mtx->entries ;
643 row = ZV_entries(rowZV) ;
644 for ( j = 0 ; j < n2 ; j++, k += inc2 ) {
645 entries[2*k] = row[2*j] ;
646 entries[2*k+1] = row[2*j+1] ;
647 }
648 return ; }
649
650 /*--------------------------------------------------------------------*/
651 /*
652 -----------------------
653 set mtx(*,jcol) = y[*]
654
655 created -- 98may01, cca
656 -----------------------
657 */
658 void
A2_setColumnDV(A2 * mtx,DV * colDV,int jcol)659 A2_setColumnDV (
660 A2 *mtx,
661 DV *colDV,
662 int jcol
663 ) {
664 double *col, *entries ;
665 int inc1, i, k, n1 ;
666 /*
667 ---------------
668 check the input
669 ---------------
670 */
671 if ( mtx == NULL || colDV == NULL || DV_size(colDV) != (n1 = mtx->n1)
672 || jcol < 0 || jcol >= mtx->n2 ) {
673 fprintf(stderr, "\n fatal error in A2_setColumnDV(%p,%p,%d)"
674 "\n bad input\n", mtx, colDV, jcol) ;
675 exit(-1) ;
676 }
677 if ( ! A2_IS_REAL(mtx) ) {
678 fprintf(stderr, "\n fatal error in A2_setColumnDV(%p,%p,%d)"
679 "\n bad type %d, must be SPOOLES_REAL\n",
680 mtx, colDV, jcol, mtx->type) ;
681 exit(-1) ;
682 }
683 k = jcol * mtx->inc2 ;
684 inc1 = mtx->inc1 ;
685 entries = mtx->entries ;
686 col = DV_entries(colDV) ;
687 for ( i = 0 ; i < n1 ; i++, k += inc1 ) {
688 entries[k] = col[i] ;
689 }
690 return ; }
691
692 /*--------------------------------------------------------------------*/
693 /*
694 -----------------------
695 set mtx(*,jcol) = y[*]
696
697 created -- 98may01, cca
698 -----------------------
699 */
700 void
A2_setColumnZV(A2 * mtx,ZV * colZV,int jcol)701 A2_setColumnZV (
702 A2 *mtx,
703 ZV *colZV,
704 int jcol
705 ) {
706 double *col, *entries ;
707 int inc1, i, k, n1 ;
708 /*
709 ---------------
710 check the input
711 ---------------
712 */
713 if ( mtx == NULL || colZV == NULL || ZV_size(colZV) != (n1 = mtx->n1)
714 || jcol < 0 || jcol >= mtx->n2 ) {
715 fprintf(stderr, "\n fatal error in A2_setColumnZV(%p,%p,%d)"
716 "\n bad input\n", mtx, colZV, jcol) ;
717 exit(-1) ;
718 }
719 if ( ! A2_IS_COMPLEX(mtx) ) {
720 fprintf(stderr, "\n fatal error in A2_setColumnZV(%p,%p,%d)"
721 "\n bad type %d, must be SPOOLES_COMPLEX\n",
722 mtx, colZV, jcol, mtx->type) ;
723 exit(-1) ;
724 }
725 k = jcol * mtx->inc2 ;
726 inc1 = mtx->inc1 ;
727 entries = mtx->entries ;
728 col = ZV_entries(colZV) ;
729 for ( i = 0 ; i < n1 ; i++, k += inc1 ) {
730 entries[2*k] = col[2*i] ;
731 entries[2*k+1] = col[2*i+1] ;
732 }
733 return ; }
734
735 /*--------------------------------------------------------------------*/
736 /*
737 -------------------------------------------------------------
738 fill the matrix with uniform random numbers in [lower, upper]
739
740 created -- 98may01, cca
741 -------------------------------------------------------------
742 */
743 void
A2_fillRandomUniform(A2 * a,double lower,double upper,int seed)744 A2_fillRandomUniform (
745 A2 *a,
746 double lower,
747 double upper,
748 int seed
749 ) {
750 double *entries ;
751 int i, inc1, inc2, j, loc, n1, n2 ;
752 Drand drand ;
753 /*
754 ---------------
755 check the input
756 ---------------
757 */
758 if ( a == NULL
759 || (n1 = a->n1) <= 0
760 || (n2 = a->n2) <= 0
761 || (inc1 = a->inc1) <= 0
762 || (inc2 = a->inc2) <= 0
763 || (entries = a->entries) == NULL ) {
764 fprintf(stderr,
765 "\n fatal error in A2_fillRandomUniform(%p,%f,%f,%d)"
766 "\n bad input\n",
767 a, lower, upper, seed) ;
768 exit(-1) ;
769 }
770 if ( ! (A2_IS_REAL(a) || A2_IS_COMPLEX(a)) ) {
771 fprintf(stderr, "\n fatal error in A2_fillRandomUniform(%p,%f,%f,%d)"
772 "\n bad type %d, must be SPOOLES_REAL or SPOOLES_COMPLEX\n",
773 a, lower, upper, seed, a->type) ;
774 exit(-1) ;
775 }
776 /*
777 ----------------
778 fill the entries
779 ----------------
780 */
781 Drand_setDefaultFields(&drand) ;
782 Drand_init(&drand) ;
783 Drand_setSeed(&drand, seed) ;
784 Drand_setUniform(&drand, lower, upper) ;
785 for ( j = 0 ; j < n2 ; j++ ) {
786 for ( i = 0 ; i < n1 ; i++ ) {
787 loc = i*inc1 + j*inc2 ;
788 if ( A2_IS_REAL(a) ) {
789 entries[loc] = Drand_value(&drand) ;
790 } else if ( A2_IS_COMPLEX(a) ) {
791 entries[2*loc] = Drand_value(&drand) ;
792 entries[2*loc+1] = Drand_value(&drand) ;
793 }
794 }
795 }
796 return ; }
797
798 /*--------------------------------------------------------------------*/
799 /*
800 -----------------------------------------------
801 fill the matrix with normal(0,1) random numbers
802
803 created -- 98may01, cca
804 -----------------------------------------------
805 */
806 void
A2_fillRandomNormal(A2 * a,double mean,double variance,int seed)807 A2_fillRandomNormal (
808 A2 *a,
809 double mean,
810 double variance,
811 int seed
812 ) {
813 double *entries ;
814 int i, inc1, inc2, j, loc, n1, n2 ;
815 Drand drand ;
816 /*
817 ---------------
818 check the input
819 ---------------
820 */
821 if ( a == NULL
822 || (n1 = a->n1) <= 0
823 || (n2 = a->n2) <= 0
824 || (inc1 = a->inc1) <= 0
825 || (inc2 = a->inc2) <= 0
826 || (entries = a->entries) == NULL ) {
827 fprintf(stderr, "\n fatal error in A2_fillRandomNormal(%p,%d)"
828 "\n bad input\n",
829 a, seed) ;
830 exit(-1) ;
831 }
832 if ( ! (A2_IS_REAL(a) || A2_IS_COMPLEX(a)) ) {
833 fprintf(stderr, "\n fatal error in A2_fillRandomNormal(%p,%f,%f,%d)"
834 "\n bad type %d, must be SPOOLES_REAL or SPOOLES_COMPLEX\n",
835 a, mean, variance, seed, a->type) ;
836 exit(-1) ;
837 }
838 /*
839 ----------------
840 fill the entries
841 ----------------
842 */
843 Drand_setDefaultFields(&drand) ;
844 Drand_init(&drand) ;
845 Drand_setSeed(&drand, seed) ;
846 Drand_setUniform(&drand, mean, variance) ;
847 for ( j = 0 ; j < n2 ; j++ ) {
848 for ( i = 0 ; i < n1 ; i++ ) {
849 loc = i*inc1 + j*inc2 ;
850 if ( A2_IS_REAL(a) ) {
851 entries[loc] = Drand_value(&drand) ;
852 } else if ( A2_IS_COMPLEX(a) ) {
853 entries[2*loc] = Drand_value(&drand) ;
854 entries[2*loc+1] = Drand_value(&drand) ;
855 }
856 }
857 }
858
859 return ; }
860
861 /*--------------------------------------------------------------------*/
862 /*
863 ----------------------------------------
864 fill the matrix with the identity matrix
865
866 created -- 98may01, cca
867 ----------------------------------------
868 */
869 void
A2_fillWithIdentity(A2 * a)870 A2_fillWithIdentity (
871 A2 *a
872 ) {
873 double *entries ;
874 int ii, inc, inc1, inc2, j, n ;
875 /*
876 ---------------
877 check the input
878 ---------------
879 */
880 if ( a == NULL
881 || (n = a->n1) <= 0
882 || n != a->n2
883 || (inc1 = a->inc1) <= 0
884 || (inc2 = a->inc2) <= 0
885 || (inc1 != 1 && inc2 != 1)
886 || (entries = a->entries) == NULL ) {
887 fprintf(stderr, "\n fatal error in A2_fillWithIdentity(%p)"
888 "\n bad input\n", a) ;
889 exit(-1) ;
890 }
891 if ( ! (A2_IS_REAL(a) || A2_IS_COMPLEX(a)) ) {
892 fprintf(stderr, "\n fatal error in A2_fillWithIdentity(%p)"
893 "\n bad type %d, must be SPOOLES_REAL or SPOOLES_COMPLEX\n",
894 a, a->type) ;
895 exit(-1) ;
896 }
897 inc = (inc1 == 1) ? inc2 : inc1 ;
898 A2_zero(a) ;
899 for ( j = 0, ii = 0 ; j < n ; j++, ii += inc + 1 ) {
900 if ( A2_IS_REAL(a) ) {
901 entries[ii] = 1.0 ;
902 } else if ( A2_IS_COMPLEX(a) ) {
903 entries[2*ii] = 1.0 ;
904 }
905 }
906 return ; }
907
908 /*--------------------------------------------------------------------*/
909 /*
910 --------------------------
911 fill the matrix with zeros
912
913 created -- 98may01, cca
914 --------------------------
915 */
916 void
A2_zero(A2 * a)917 A2_zero (
918 A2 *a
919 ) {
920 double *entries ;
921 int i, inc1, inc2, j, loc, n1, n2 ;
922 /*
923 ---------------
924 check the input
925 ---------------
926 */
927 if ( a == NULL
928 || (n1 = a->n1) <= 0
929 || (n2 = a->n2) <= 0
930 || (inc1 = a->inc1) <= 0
931 || (inc2 = a->inc2) <= 0
932 || (entries = a->entries) == NULL ) {
933 fprintf(stderr, "\n fatal error in A2_zero(%p)"
934 "\n bad input\n", a) ;
935 exit(-1) ;
936 }
937 if ( ! (A2_IS_REAL(a) || A2_IS_COMPLEX(a)) ) {
938 fprintf(stderr, "\n fatal error in A2_zero(%p)"
939 "\n bad type %d, must be SPOOLES_REAL or SPOOLES_COMPLEX\n",
940 a, a->type) ;
941 exit(-1) ;
942 }
943 for ( j = 0 ; j < n2 ; j++ ) {
944 for ( i = 0 ; i < n1 ; i++ ) {
945 loc =i*inc1 + j*inc2 ;
946 if ( A2_IS_REAL(a) ) {
947 entries[loc] = 0.0 ;
948 } else if ( A2_IS_COMPLEX(a) ) {
949 entries[2*loc] = 0.0 ;
950 entries[2*loc+1] = 0.0 ;
951 }
952 }
953 }
954
955 return ; }
956
957 /*--------------------------------------------------------------------*/
958 /*
959 ----------------------------
960 copy one matrix into another
961 A := B
962
963 created -- 98may01, cca
964 ----------------------------
965 */
966 void
A2_copy(A2 * A,A2 * B)967 A2_copy (
968 A2 *A,
969 A2 *B
970 ) {
971 double *entA, *entB ;
972 int inc1A, inc1B, inc2A, inc2B, irow, jcol, locA, locB,
973 ncol, ncolA, ncolB, nrow, nrowA, nrowB ;
974 /*
975 ---------------
976 check the input
977 ---------------
978 */
979 if ( A == NULL
980 || (nrowA = A->n1) < 0
981 || (ncolA = A->n2) < 0
982 || (inc1A = A->inc1) <= 0
983 || (inc2A = A->inc2) <= 0
984 || (entA = A->entries) == NULL
985 || B == NULL
986 || (nrowB = B->n1) < 0
987 || (ncolB = B->n2) < 0
988 || (inc1B = B->inc1) <= 0
989 || (inc2B = B->inc2) <= 0
990 || (entB = B->entries) == NULL ) {
991 fprintf(stderr, "\n fatal error in A2_copy(%p,%p)"
992 "\n bad input\n", A, B) ;
993 if ( A != NULL ) {
994 fprintf(stderr, "\n\n first A2 object") ;
995 A2_writeStats(A, stderr) ;
996 }
997 if ( B != NULL ) {
998 fprintf(stderr, "\n\n second A2 object") ;
999 A2_writeStats(B, stderr) ;
1000 }
1001 exit(-1) ;
1002 }
1003 if ( ! (A2_IS_REAL(A) || A2_IS_COMPLEX(A)) ) {
1004 fprintf(stderr, "\n fatal error in A2_copy(%p,%p)"
1005 "\n bad type %d, must be SPOOLES_REAL or SPOOLES_COMPLEX\n",
1006 A, B, A->type) ;
1007 exit(-1) ;
1008 }
1009 if ( ! (A2_IS_REAL(B) || A2_IS_COMPLEX(B)) ) {
1010 fprintf(stderr, "\n fatal error in A2_copy(%p,%p)"
1011 "\n bad type %d, must be SPOOLES_REAL or SPOOLES_COMPLEX\n",
1012 A, B, B->type) ;
1013 exit(-1) ;
1014 }
1015 if ( A->type != B->type ) {
1016 fprintf(stderr, "\n fatal error in A2_copy(%p,%p)"
1017 "\n A's type %d, B's type = %d, must be the same\n",
1018 A, B, A->type, B->type) ;
1019 exit(-1) ;
1020 }
1021 nrow = (nrowA <= nrowB) ? nrowA : nrowB ;
1022 ncol = (ncolA <= ncolB) ? ncolA : ncolB ;
1023 if ( A2_IS_REAL(A) ) {
1024 if ( inc1A == 1 && inc1B == 1 ) {
1025 double *colA = entA, *colB = entB ;
1026 for ( jcol = 0 ; jcol < ncol ; jcol++ ) {
1027 for ( irow = 0 ; irow < nrow ; irow++ ) {
1028 colA[irow] = colB[irow] ;
1029 }
1030 colA += inc2A ;
1031 colB += inc2B ;
1032 }
1033 } else if ( inc2A == 1 && inc2B == 1 ) {
1034 double *rowA = entA, *rowB = entB ;
1035 for ( irow = 0 ; irow < nrow ; irow++ ) {
1036 for ( jcol = 0 ; jcol < ncol ; jcol++ ) {
1037 rowA[jcol] = rowB[jcol] ;
1038 }
1039 rowA += 2*inc1A ;
1040 }
1041 } else {
1042 for ( irow = 0 ; irow < nrow ; irow++ ) {
1043 for ( jcol = 0 ; jcol < ncol ; jcol++ ) {
1044 locA = irow*inc1A + jcol*inc2A ;
1045 locB = irow*inc1B + jcol*inc2B ;
1046 entA[locA] = entB[locB] ;
1047 }
1048 }
1049 }
1050 } else if ( A2_IS_COMPLEX(A) ) {
1051 if ( inc1A == 1 && inc1B == 1 ) {
1052 double *colA = entA, *colB = entB ;
1053 for ( jcol = 0 ; jcol < ncol ; jcol++ ) {
1054 for ( irow = 0 ; irow < nrow ; irow++ ) {
1055 colA[2*irow] = colB[2*irow] ;
1056 colA[2*irow+1] = colB[2*irow+1] ;
1057 }
1058 colA += 2*inc2A ;
1059 colB += 2*inc2B ;
1060 }
1061 } else if ( inc2A == 1 && inc2B == 1 ) {
1062 double *rowA = entA, *rowB = entB ;
1063 for ( irow = 0 ; irow < nrow ; irow++ ) {
1064 for ( jcol = 0 ; jcol < ncol ; jcol++ ) {
1065 rowA[2*jcol] = rowB[2*jcol] ;
1066 rowA[2*jcol+1] = rowB[2*jcol+1] ;
1067 }
1068 rowA += 2*inc1A ;
1069 rowB += 2*inc1B ;
1070 }
1071 } else {
1072 for ( irow = 0 ; irow < nrow ; irow++ ) {
1073 for ( jcol = 0 ; jcol < ncol ; jcol++ ) {
1074 locA = irow*inc1A + jcol*inc2A ;
1075 locB = irow*inc1B + jcol*inc2B ;
1076 entA[2*locA] = entB[2*locB] ;
1077 entA[2*locA+1] = entB[2*locB+1] ;
1078 }
1079 }
1080 }
1081 }
1082 return ; }
1083
1084 /*--------------------------------------------------------------------*/
1085 /*
1086 --------------------------------
1087 subtract one matrix from another
1088
1089 A := A - B
1090
1091 created -- 98may01, cca
1092 ----------------------------
1093 */
1094 void
A2_sub(A2 * A,A2 * B)1095 A2_sub (
1096 A2 *A,
1097 A2 *B
1098 ) {
1099 double *entA, *entB ;
1100 int inc1A, inc1B, inc2A, inc2B, irow, jcol, locA, locB,
1101 ncol, ncolA, ncolB, nrow, nrowA, nrowB ;
1102 /*
1103 ---------------
1104 check the input
1105 ---------------
1106 */
1107 if ( A == NULL
1108 || B == NULL
1109 || (nrowA = A->n1) <= 0
1110 || (ncolA = A->n2) <= 0
1111 || (inc1A = A->inc1) <= 0
1112 || (inc2A = A->inc2) <= 0
1113 || (nrowB = B->n1) <= 0
1114 || (ncolB = B->n2) <= 0
1115 || (inc1B = B->inc1) <= 0
1116 || (inc2B = B->inc2) <= 0
1117 || (entA = A->entries) == NULL
1118 || (entB = B->entries) == NULL ) {
1119 fprintf(stderr, "\n fatal error in A2_sub(%p,%p)"
1120 "\n bad input\n", A, B) ;
1121 if ( A != NULL ) {
1122 fprintf(stderr, "\n\n first A2 object") ;
1123 A2_writeStats(A, stderr) ;
1124 }
1125 if ( B != NULL ) {
1126 fprintf(stderr, "\n\n second A2 object") ;
1127 A2_writeStats(B, stderr) ;
1128 }
1129 exit(-1) ;
1130 }
1131 if ( ! (A2_IS_REAL(A) || A2_IS_COMPLEX(A)) ) {
1132 fprintf(stderr, "\n fatal error in A2_sub(%p,%p)"
1133 "\n bad type %d, must be SPOOLES_REAL or SPOOLES_COMPLEX\n",
1134 A, B, A->type) ;
1135 exit(-1) ;
1136 }
1137 if ( ! (A2_IS_REAL(B) || A2_IS_COMPLEX(B)) ) {
1138 fprintf(stderr, "\n fatal error in A2_sub(%p,%p)"
1139 "\n bad type %d, must be SPOOLES_REAL or SPOOLES_COMPLEX\n",
1140 A, B, B->type) ;
1141 exit(-1) ;
1142 }
1143 if ( A->type != B->type ) {
1144 fprintf(stderr, "\n fatal error in A2_sub(%p,%p)"
1145 "\n A's type %d, B's type = %d, must be the same\n",
1146 A, B, A->type, B->type) ;
1147 exit(-1) ;
1148 }
1149 /*
1150 fprintf(stdout, "\n debug : A") ;
1151 A2_writeForHumanEye(A, stdout) ;
1152 fprintf(stdout, "\n debug : B") ;
1153 A2_writeForHumanEye(B, stdout) ;
1154 */
1155 nrow = (nrowA <= nrowB) ? nrowA : nrowB ;
1156 ncol = (ncolA <= ncolB) ? ncolA : ncolB ;
1157 if ( A2_IS_REAL(A) ) {
1158 for ( irow = 0 ; irow < nrow ; irow++ ) {
1159 for ( jcol = 0 ; jcol < ncol ; jcol++ ) {
1160 locA = irow*inc1A + jcol*inc2A ;
1161 locB = irow*inc1B + jcol*inc2B ;
1162 entA[locA] -= entB[locB] ;
1163 }
1164 }
1165 } else if ( A2_IS_COMPLEX(A) ) {
1166 for ( irow = 0 ; irow < nrow ; irow++ ) {
1167 for ( jcol = 0 ; jcol < ncol ; jcol++ ) {
1168 locA = irow*inc1A + jcol*inc2A ;
1169 locB = irow*inc1B + jcol*inc2B ;
1170 entA[2*locA] -= entB[2*locB] ;
1171 entA[2*locA+1] -= entB[2*locB+1] ;
1172 }
1173 }
1174 }
1175 return ; }
1176
1177 /*--------------------------------------------------------------------*/
1178 /*
1179 ---------------------------
1180 swap two rows of the matrix
1181
1182 created -- 98may01, cca
1183 ---------------------------
1184 */
1185 void
A2_swapRows(A2 * a,int irow1,int irow2)1186 A2_swapRows (
1187 A2 *a,
1188 int irow1,
1189 int irow2
1190 ) {
1191 double temp ;
1192 double *row1, *row2 ;
1193 int inc2, j, k, n2 ;
1194 /*
1195 -----------
1196 check input
1197 -----------
1198 */
1199 if ( a == NULL
1200 || irow1 < 0 || irow1 >= a->n1
1201 || irow2 < 0 || irow2 >= a->n1 ) {
1202 fprintf(stderr,
1203 "\n fatal error in A2_swapRows(%p,%d,%d)"
1204 "\n bad input\n", a, irow1, irow2) ;
1205 exit(-1) ;
1206 }
1207 if ( a->n1 <= 0
1208 || a->inc1 <= 0
1209 || (n2 = a->n2) <= 0
1210 || (inc2 = a->inc2) <= 0
1211 || a->entries == NULL ) {
1212 fprintf(stderr,
1213 "\n fatal error in A2_swapRows(%p,%d,%d)"
1214 "\n bad structure\n", a, irow1, irow2) ;
1215 exit(-1) ;
1216 }
1217 if ( ! (A2_IS_REAL(a) || A2_IS_COMPLEX(a)) ) {
1218 fprintf(stderr, "\n fatal error in A2_swapRows(%p,%d,%d)"
1219 "\n bad type %d, must be SPOOLES_REAL or SPOOLES_COMPLEX\n",
1220 a, irow1, irow2, a->type) ;
1221 exit(-1) ;
1222 }
1223 if ( irow1 == irow2 ) {
1224 return ;
1225 }
1226 if ( A2_IS_REAL(a) ) {
1227 row1 = a->entries + irow1*a->inc1 ;
1228 row2 = a->entries + irow2*a->inc1 ;
1229 if ( inc2 == 1 ) {
1230 for ( j = 0 ; j < n2 ; j++ ) {
1231 temp = row1[j] ;
1232 row1[j] = row2[j] ;
1233 row2[j] = temp ;
1234 }
1235 } else {
1236 for ( j = 0, k = 0 ; j < n2 ; j++, k += inc2 ) {
1237 temp = row1[k] ;
1238 row1[k] = row2[k] ;
1239 row2[k] = temp ;
1240 }
1241 }
1242 } else if ( A2_IS_COMPLEX(a) ) {
1243 row1 = a->entries + 2*irow1*a->inc1 ;
1244 row2 = a->entries + 2*irow2*a->inc1 ;
1245 if ( inc2 == 1 ) {
1246 for ( j = 0 ; j < n2 ; j++ ) {
1247 temp = row1[2*j] ;
1248 row1[2*j] = row2[2*j] ;
1249 row2[2*j] = temp ;
1250 temp = row1[2*j+1] ;
1251 row1[2*j+1] = row2[2*j+1] ;
1252 row2[2*j+1] = temp ;
1253 }
1254 } else {
1255 for ( j = 0, k = 0 ; j < n2 ; j++, k += inc2 ) {
1256 temp = row1[2*k] ;
1257 row1[2*k] = row2[2*k] ;
1258 row2[2*k] = temp ;
1259 temp = row1[2*k+1] ;
1260 row1[2*k+1] = row2[2*k+1] ;
1261 row2[2*k+1] = temp ;
1262 }
1263 }
1264 }
1265 return ; }
1266
1267 /*--------------------------------------------------------------------*/
1268 /*
1269 ------------------------------
1270 swap two columns of the matrix
1271
1272 created -- 98may01, cca
1273 ------------------------------
1274 */
1275 void
A2_swapColumns(A2 * a,int jcol1,int jcol2)1276 A2_swapColumns (
1277 A2 *a,
1278 int jcol1,
1279 int jcol2
1280 ) {
1281 double temp ;
1282 double *col1, *col2 ;
1283 int i, inc1, k, n1 ;
1284 /*
1285 -----------
1286 check input
1287 -----------
1288 */
1289 if ( a == NULL
1290 || jcol1 < 0 || jcol1 >= a->n2
1291 || jcol2 < 0 || jcol2 >= a->n2 ) {
1292 fprintf(stderr,
1293 "\n fatal error in A2_swapColumns(%p,%d,%d)"
1294 "\n bad input\n", a, jcol1, jcol2) ;
1295 exit(-1) ;
1296 }
1297 if ( (n1 = a->n1) <= 0
1298 || (inc1 = a->inc1) <= 0
1299 || a->n2 <= 0
1300 || a->inc2 <= 0
1301 || a->entries == NULL ) {
1302 fprintf(stderr,
1303 "\n fatal error in A2_swapColumns(%p,%d,%d)"
1304 "\n bad structure\n", a, jcol1, jcol2) ;
1305 exit(-1) ;
1306 }
1307 if ( ! (A2_IS_REAL(a) || A2_IS_COMPLEX(a)) ) {
1308 fprintf(stderr, "\n fatal error in A2_swapColumns(%p,%d,%d)"
1309 "\n bad type %d, must be SPOOLES_REAL or SPOOLES_COMPLEX\n",
1310 a, jcol1, jcol2, a->type) ;
1311 exit(-1) ;
1312 }
1313 if ( jcol1 == jcol2 ) {
1314 return ;
1315 }
1316 if ( A2_IS_REAL(a) ) {
1317 col1 = a->entries + jcol1*a->inc2 ;
1318 col2 = a->entries + jcol2*a->inc2 ;
1319 if ( inc1 == 1 ) {
1320 for ( i = 0 ; i < n1 ; i++ ) {
1321 temp = col1[i] ;
1322 col1[i] = col2[i] ;
1323 col2[i] = temp ;
1324 }
1325 } else {
1326 for ( i = 0, k = 0 ; i < n1 ; i++, k += inc1 ) {
1327 temp = col1[k] ;
1328 col1[k] = col2[k] ;
1329 col2[k] = temp ;
1330 }
1331 }
1332 } else if ( A2_IS_COMPLEX(a) ) {
1333 col1 = a->entries + 2*jcol1*a->inc2 ;
1334 col2 = a->entries + 2*jcol2*a->inc2 ;
1335 if ( inc1 == 1 ) {
1336 for ( i = 0 ; i < n1 ; i++ ) {
1337 temp = col1[2*i] ;
1338 col1[2*i] = col2[2*i] ;
1339 col2[2*i] = temp ;
1340 temp = col1[2*i+1] ;
1341 col1[2*i+1] = col2[2*i+1] ;
1342 col2[2*i+1] = temp ;
1343 }
1344 } else {
1345 for ( i = 0, k = 0 ; i < n1 ; i++, k += inc1 ) {
1346 temp = col1[2*k] ;
1347 col1[2*k] = col2[2*k] ;
1348 col2[2*k] = temp ;
1349 temp = col1[2*k+1] ;
1350 col1[2*k+1] = col2[2*k+1] ;
1351 col2[2*k+1] = temp ;
1352 }
1353 }
1354 }
1355 return ; }
1356
1357 /*--------------------------------------------------------------------*/
1358