1 /* norms.c */
2
3 #include "../A2.h"
4
5 #define CAUTIOUS 1
6
7 /*--------------------------------------------------------------------*/
8 /*
9 -------------------------------------
10 return the entry of maximum magnitude
11
12 created -- 98apr15, cca
13 -------------------------------------
14 */
15 double
A2_maxabs(A2 * a)16 A2_maxabs (
17 A2 *a
18 ) {
19 double maxval, val ;
20 double *entries, *row ;
21 int inc1, inc2, irow, jcol, kk, n1, n2 ;
22 /*
23 ---------------
24 check the input
25 ---------------
26 */
27 if ( a == NULL
28 || (n1 = a->n1) < 0
29 || (n2 = a->n2) < 0
30 || (inc1 = a->inc1) < 0
31 || (inc2 = a->inc2) < 0 ) {
32 fprintf(stderr, "\n fatal error in A2_maxabs(%p)"
33 "\n bad input\n", a ) ;
34 exit(-1) ;
35 }
36 if ( ! (A2_IS_REAL(a) || A2_IS_COMPLEX(a)) ) {
37 fprintf(stderr, "\n fatal error in A2_maxabs(%p)"
38 "\n bad type %d, must be SPOOLES_REAL or SPOOLES_COMPLEX\n",
39 a, a->type) ;
40 exit(-1) ;
41 }
42 entries = a->entries ;
43 maxval = 0.0 ;
44 row = entries ;
45 if ( A2_IS_REAL(a) ) {
46 for ( irow = 0 ; irow < n1 ; irow++ ) {
47 for ( jcol = 0, kk = 0 ; jcol < n2 ; jcol++, kk += inc2 ) {
48 val = fabs(row[kk]) ;
49 if ( maxval < val ) {
50 maxval = val ;
51 }
52 }
53 row += inc1 ;
54 }
55 } else if ( A2_IS_COMPLEX(a) ) {
56 for ( irow = 0 ; irow < n1 ; irow++ ) {
57 for ( jcol = 0, kk = 0 ; jcol < n2 ; jcol++, kk += inc2 ) {
58 val = Zabs(row[2*kk], row[2*kk+1]) ;
59 if ( maxval < val ) {
60 maxval = val ;
61 }
62 }
63 row += inc1 ;
64 }
65 }
66 return(maxval) ; }
67
68 /*--------------------------------------------------------------------*/
69 /*
70 ---------------------------------------
71 return the frobenius norm of the matrix
72
73 created -- 98apr15, cca
74 ---------------------------------------
75 */
76 double
A2_frobNorm(A2 * mtx)77 A2_frobNorm (
78 A2 *mtx
79 ) {
80 double norm ;
81 int ncol, nrow ;
82 /*
83 ---------------
84 check the input
85 ---------------
86 */
87 if ( mtx == NULL ) {
88 fprintf(stderr,
89 "\n fatal error in A2_frobNorm(%p) "
90 "\n bad input\n", mtx) ;
91 exit(-1) ;
92 }
93 if ( ! (A2_IS_REAL(mtx) || A2_IS_COMPLEX(mtx)) ) {
94 fprintf(stderr, "\n fatal error in A2_frobNorm(%p)"
95 "\n bad type %d, must be SPOOLES_REAL or SPOOLES_COMPLEX\n",
96 mtx, mtx->type) ;
97 exit(-1) ;
98 }
99 if ( (nrow = mtx->n1) <= 0 || (ncol = mtx->n2) <= 0 ) {
100 return(0.0) ;
101 }
102 norm = 0 ;
103 if ( A2_IS_REAL(mtx) ) {
104 if ( mtx->inc1 == 1 ) {
105 double *col ;
106 int inc2 = mtx->inc2, irow, jcol ;
107
108 for ( jcol = 0, col = mtx->entries ;
109 jcol < ncol ;
110 jcol++, col += inc2 ) {
111 for ( irow = 0 ; irow < nrow ; irow++ ) {
112 norm += col[irow]*col[irow] ;
113 }
114 }
115 } else if ( mtx->inc2 == 1 ) {
116 double *row ;
117 int inc1 = mtx->inc1, irow, jcol ;
118
119 for ( irow = 0, row = mtx->entries ;
120 irow < nrow ;
121 irow++, row += inc1 ) {
122 for ( jcol = 0 ; jcol < ncol ; jcol++ ) {
123 norm += row[jcol]*row[jcol] ;
124 }
125 }
126 } else {
127 double *entries = mtx->entries ;
128 int inc1 = mtx->inc1, inc2 = mtx->inc2, irow, jcol, loc ;
129
130 for ( irow = 0 ; irow < nrow ; irow++ ) {
131 for ( jcol = 0 ; jcol < ncol ; jcol++ ) {
132 loc = irow*inc1 + jcol*inc2 ;
133 norm += entries[loc]*entries[loc] ;
134 }
135 }
136 }
137 } else if ( A2_IS_COMPLEX(mtx) ) {
138 if ( mtx->inc1 == 1 ) {
139 double *col ;
140 int inc2 = mtx->inc2, irow, jcol ;
141
142 for ( jcol = 0, col = mtx->entries ;
143 jcol < ncol ;
144 jcol++, col += 2*inc2 ) {
145 for ( irow = 0 ; irow < nrow ; irow++ ) {
146 norm += col[2*irow]*col[2*irow]
147 + col[2*irow+1]*col[2*irow+1] ;
148 }
149 }
150 } else if ( mtx->inc2 == 1 ) {
151 double *row ;
152 int inc1 = mtx->inc1, irow, jcol ;
153
154 for ( irow = 0, row = mtx->entries ;
155 irow < nrow ;
156 irow++, row += 2*inc1 ) {
157 for ( jcol = 0 ; jcol < ncol ; jcol++ ) {
158 norm += row[2*jcol]*row[2*jcol]
159 + row[2*jcol+1]*row[2*jcol+1] ;
160 }
161 }
162 } else {
163 double *entries = mtx->entries ;
164 int inc1 = mtx->inc1, inc2 = mtx->inc2, irow, jcol, loc ;
165
166 for ( irow = 0 ; irow < nrow ; irow++ ) {
167 for ( jcol = 0 ; jcol < ncol ; jcol++ ) {
168 loc = irow*inc1 + jcol*inc2 ;
169 norm += entries[2*loc]*entries[2*loc]
170 + entries[2*loc+1]*entries[2*loc+1] ;
171 }
172 }
173 }
174 }
175 norm = sqrt(norm) ;
176
177 return(norm) ; }
178
179 /*--------------------------------------------------------------------*/
180 /*
181 ---------------------------------
182 return the one-norm of the matrix
183
184 created -- 98apr15, cca
185 ---------------------------------
186 */
187 double
A2_oneNorm(A2 * mtx)188 A2_oneNorm (
189 A2 *mtx
190 ) {
191 double norm ;
192 int ncol, nrow ;
193 /*
194 ---------------
195 check the input
196 ---------------
197 */
198 if ( mtx == NULL ) {
199 fprintf(stderr, "\n fatal error in A2_oneNorm(%p) "
200 "\n bad input\n", mtx) ;
201 exit(-1) ;
202 }
203 if ( ! (A2_IS_REAL(mtx) || A2_IS_COMPLEX(mtx)) ) {
204 fprintf(stderr, "\n fatal error in A2_oneNorm(%p)"
205 "\n bad type %d, must be SPOOLES_REAL or SPOOLES_COMPLEX\n",
206 mtx, mtx->type) ;
207 exit(-1) ;
208 }
209 if ( (nrow = mtx->n1) <= 0 || (ncol = mtx->n2) <= 0 ) {
210 return(0.0) ;
211 }
212 norm = 0.0 ;
213 if ( A2_IS_REAL(mtx) ) {
214 if ( mtx->inc1 == 1 ) {
215 double sum ;
216 double *col ;
217 int inc2 = mtx->inc2, irow, jcol ;
218
219 for ( jcol = 0, col = mtx->entries ;
220 jcol < ncol ;
221 jcol++, col += inc2 ) {
222 sum = 0.0 ;
223 for ( irow = 0 ; irow < nrow ; irow++ ) {
224 sum += fabs(col[irow]) ;
225 }
226 if ( norm < sum ) {
227 norm = sum ;
228 }
229 }
230 } else {
231 double sum ;
232 double *col ;
233 int inc1 = mtx->inc1, inc2 = mtx->inc2, irow, jcol, kk ;
234
235 for ( jcol = 0, col = mtx->entries ;
236 jcol < ncol ;
237 jcol++, col += inc2 ) {
238 sum = 0.0 ;
239 for ( irow = 0, kk = 0 ; irow < nrow ; irow++, kk += inc1 ) {
240 sum += fabs(col[kk]) ;
241 }
242 if ( norm < sum ) {
243 norm = sum ;
244 }
245 }
246 }
247 } else if ( A2_IS_COMPLEX(mtx) ) {
248 if ( mtx->inc1 == 1 ) {
249 double sum ;
250 double *col ;
251 int inc2 = mtx->inc2, irow, jcol ;
252
253 for ( jcol = 0, col = mtx->entries ;
254 jcol < ncol ;
255 jcol++, col += 2*inc2 ) {
256 sum = 0.0 ;
257 for ( irow = 0 ; irow < nrow ; irow++ ) {
258 sum += Zabs(col[2*irow], col[2*irow+1]) ;
259 }
260 if ( norm < sum ) {
261 norm = sum ;
262 }
263 }
264 } else {
265 double sum ;
266 double *col ;
267 int inc1 = mtx->inc1, inc2 = mtx->inc2, irow, jcol, kk ;
268
269 for ( jcol = 0, col = mtx->entries ;
270 jcol < ncol ;
271 jcol++, col += 2*inc2 ) {
272 sum = 0.0 ;
273 for ( irow = 0, kk = 0 ; irow < nrow ; irow++, kk += inc1 ) {
274 sum += Zabs(col[2*kk], col[2*kk+1]) ;
275 }
276 if ( norm < sum ) {
277 norm = sum ;
278 }
279 }
280 }
281 }
282 return(norm) ; }
283
284 /*--------------------------------------------------------------------*/
285 /*
286 --------------------------------------
287 return the infinity-norm of the matrix
288
289 created -- 98apr15, cca
290 --------------------------------------
291 */
292 double
A2_infinityNorm(A2 * mtx)293 A2_infinityNorm (
294 A2 *mtx
295 ) {
296 double norm ;
297 int ncol, nrow ;
298 /*
299 ---------------
300 check the input
301 ---------------
302 */
303 if ( mtx == NULL ) {
304 fprintf(stderr, "\n fatal error in A2_infinityNorm(%p) "
305 "\n bad input\n", mtx) ;
306 exit(-1) ;
307 }
308 if ( ! (A2_IS_REAL(mtx) || A2_IS_COMPLEX(mtx)) ) {
309 fprintf(stderr, "\n fatal error in A2_infinityNorm(%p)"
310 "\n bad type %d, must be SPOOLES_REAL or SPOOLES_COMPLEX\n",
311 mtx, mtx->type) ;
312 exit(-1) ;
313 }
314 if ( (nrow = mtx->n1) <= 0 || (ncol = mtx->n2) <= 0 ) {
315 return(0.0) ;
316 }
317 norm = 0.0 ;
318 if ( A2_IS_REAL(mtx) ) {
319 if ( mtx->inc2 == 1 ) {
320 double sum ;
321 double *row = mtx->entries ;
322 int inc1 = mtx->inc1, irow, jcol ;
323
324 for ( irow = 0 ; irow < nrow ; irow++, row += inc1 ) {
325 for ( jcol = 0, sum = 0.0 ; jcol < ncol ; jcol++ ) {
326 sum += fabs(row[jcol]) ;
327 }
328 if ( norm < sum ) {
329 norm = sum ;
330 }
331 }
332 } else {
333 double sum ;
334 double *row = mtx->entries ;
335 int inc1 = mtx->inc1, inc2 = mtx->inc2, irow, jcol, kk ;
336
337 for ( irow = 0 ; irow < nrow ; irow++, row += inc1 ) {
338 sum = 0.0 ;
339 for ( jcol = 0, kk = 0 ; jcol < ncol ; jcol++, kk += inc2 ) {
340 sum += fabs(row[kk]) ;
341 }
342 if ( norm < sum ) {
343 norm = sum ;
344 }
345 }
346 }
347 } else if ( A2_IS_COMPLEX(mtx) ) {
348 if ( mtx->inc2 == 1 ) {
349 double sum ;
350 double *row ;
351 int inc1 = mtx->inc1, irow, jcol ;
352
353 for ( irow = 0, row = mtx->entries ;
354 irow < nrow ;
355 irow++, row += 2*inc1 ) {
356 sum = 0.0 ;
357 for ( jcol = 0 ; jcol < ncol ; jcol++ ) {
358 sum += Zabs(row[2*jcol], row[2*jcol+1]) ;
359 }
360 if ( norm < sum ) {
361 norm = sum ;
362 }
363 }
364 } else {
365 double sum ;
366 double *row ;
367 int inc1 = mtx->inc1, inc2 = mtx->inc2, irow, jcol, kk ;
368
369 for ( irow = 0, row = mtx->entries ;
370 irow < nrow ;
371 irow++, row += 2*inc1 ) {
372 sum = 0.0 ;
373 for ( jcol = 0, kk = 0 ; jcol < ncol ; jcol++, kk += inc2 ) {
374 sum += Zabs(row[2*kk], row[2*kk+1]) ;
375 }
376 if ( norm < sum ) {
377 norm = sum ;
378 }
379 }
380 }
381 }
382 return(norm) ; }
383
384 /*--------------------------------------------------------------------*/
385 /*
386 ----------------------------------
387 return the one-norm of column jcol
388
389 created -- 98apr15, cca
390 ----------------------------------
391 */
392 double
A2_oneNormOfColumn(A2 * mtx,int jcol)393 A2_oneNormOfColumn (
394 A2 *mtx,
395 int jcol
396 ) {
397 double sum ;
398 double *col ;
399 int inc1, irow, kk, nrow ;
400 /*
401 ---------------
402 check the input
403 ---------------
404 */
405 if ( mtx == NULL || mtx->entries == NULL
406 || jcol < 0 || jcol > mtx->n2 ) {
407 fprintf(stderr, "\n fatal error in A2_oneNormOfColumn(%p,%d)"
408 "\n bad input\n", mtx, jcol) ;
409 exit(-1) ;
410 }
411 if ( ! (A2_IS_REAL(mtx) || A2_IS_COMPLEX(mtx)) ) {
412 fprintf(stderr, "\n fatal error in A2_oneNormOfColumn(%p,%d)"
413 "\n bad type %d, must be SPOOLES_REAL or SPOOLES_COMPLEX\n",
414 mtx, jcol, mtx->type) ;
415 exit(-1) ;
416 }
417 nrow = mtx->n1 ;
418 sum = 0.0 ;
419 if ( A2_IS_REAL(mtx) ) {
420 col = mtx->entries + jcol*mtx->inc2 ;
421 if ( (inc1 = mtx->inc1) == 1 ) {
422 for ( irow = 0 ; irow < nrow ; irow++ ) {
423 sum += fabs(col[irow]) ;
424 }
425 } else {
426 for ( irow = 0, kk = 0 ; irow < nrow ; irow++, kk += inc1 ) {
427 sum += fabs(col[kk]) ;
428 }
429 }
430 } else if ( A2_IS_COMPLEX(mtx) ) {
431 col = mtx->entries + 2*jcol*mtx->inc2 ;
432 if ( (inc1 = mtx->inc1) == 1 ) {
433 for ( irow = 0 ; irow < nrow ; irow++ ) {
434 sum += Zabs(col[2*irow], col[2*irow+1]) ;
435 }
436 } else {
437 for ( irow = 0, kk = 0 ; irow < nrow ; irow++, kk += inc1 ) {
438 sum += Zabs(col[2*kk], col[2*kk+1]) ;
439 }
440 }
441 }
442 return(sum) ; }
443
444 /*--------------------------------------------------------------------*/
445 /*
446 ----------------------------------
447 return the two-norm of column jcol
448
449 created -- 98apr15, cca
450 ----------------------------------
451 */
452 double
A2_twoNormOfColumn(A2 * mtx,int jcol)453 A2_twoNormOfColumn (
454 A2 *mtx,
455 int jcol
456 ) {
457 double sum ;
458 double *col ;
459 int inc1, irow, kk, nrow ;
460 /*
461 ---------------
462 check the input
463 ---------------
464 */
465 if ( mtx == NULL || mtx->entries == NULL
466 || jcol < 0 || jcol > mtx->n2 ) {
467 fprintf(stderr, "\n fatal error in A2_twoNormOfColumn(%p,%d)"
468 "\n bad input\n", mtx, jcol) ;
469 exit(-1) ;
470 }
471 if ( ! (A2_IS_REAL(mtx) || A2_IS_COMPLEX(mtx)) ) {
472 fprintf(stderr, "\n fatal error in A2_twoNormOfColumn(%p,%d)"
473 "\n bad type %d, must be SPOOLES_REAL or SPOOLES_COMPLEX\n",
474 mtx, jcol, mtx->type) ;
475 exit(-1) ;
476 }
477 nrow = mtx->n1 ;
478 sum = 0.0 ;
479 if ( A2_IS_REAL(mtx) ) {
480 col = mtx->entries + jcol*mtx->inc2 ;
481 if ( (inc1 = mtx->inc1) == 1 ) {
482 for ( irow = 0 ; irow < nrow ; irow++ ) {
483 sum += col[irow]*col[irow] ;
484 }
485 } else {
486 for ( irow = 0, kk = 0 ; irow < nrow ; irow++, kk += inc1 ) {
487 sum += col[kk]*col[kk] ;
488 }
489 }
490 } else if ( A2_IS_COMPLEX(mtx) ) {
491 col = mtx->entries + 2*jcol*mtx->inc2 ;
492 if ( (inc1 = mtx->inc1) == 1 ) {
493 for ( irow = 0 ; irow < nrow ; irow++ ) {
494 sum += col[2*irow]*col[2*irow] + col[2*irow+1]*col[2*irow+1] ;
495 }
496 } else {
497 for ( irow = 0, kk = 0 ; irow < nrow ; irow++, kk += inc1 ) {
498 sum += col[2*kk]*col[2*kk] + col[2*kk+1]*col[2*kk+1] ;
499 }
500 }
501 }
502 sum = sqrt(sum) ;
503
504 return(sum) ; }
505
506 /*--------------------------------------------------------------------*/
507 /*
508 ---------------------------------------
509 return the infinity-norm of column jcol
510
511 created -- 98apr15, cca
512 ---------------------------------------
513 */
514 double
A2_infinityNormOfColumn(A2 * mtx,int jcol)515 A2_infinityNormOfColumn (
516 A2 *mtx,
517 int jcol
518 ) {
519 double norm, val ;
520 double *col ;
521 int inc1, irow, kk, nrow ;
522 /*
523 ---------------
524 check the input
525 ---------------
526 */
527 if ( mtx == NULL || mtx->entries == NULL
528 || jcol < 0 || jcol > mtx->n2 ) {
529 fprintf(stderr, "\n fatal error in A2_infinityNormOfColumn(%p,%d)"
530 "\n bad input\n", mtx, jcol) ;
531 exit(-1) ;
532 }
533 if ( ! (A2_IS_REAL(mtx) || A2_IS_COMPLEX(mtx)) ) {
534 fprintf(stderr, "\n fatal error in A2_infinityNormOfColumn(%p,%d)"
535 "\n bad type %d, must be SPOOLES_REAL or SPOOLES_COMPLEX\n",
536 mtx, jcol, mtx->type) ;
537 exit(-1) ;
538 }
539 nrow = mtx->n1 ;
540 norm = 0.0 ;
541 if ( A2_IS_REAL(mtx) ) {
542 col = mtx->entries + jcol*mtx->inc2 ;
543 if ( (inc1 = mtx->inc1) == 1 ) {
544 for ( irow = 0 ; irow < nrow ; irow++ ) {
545 val = fabs(col[irow]) ;
546 if ( norm < val ) {
547 norm = val ;
548 }
549 }
550 } else {
551 for ( irow = 0, kk = 0 ; irow < nrow ; irow++, kk += inc1 ) {
552 val = fabs(col[kk]) ;
553 if ( norm < val ) {
554 norm = val ;
555 }
556 }
557 }
558 } else if ( A2_IS_COMPLEX(mtx) ) {
559 col = mtx->entries + 2*jcol*mtx->inc2 ;
560 if ( (inc1 = mtx->inc1) == 1 ) {
561 for ( irow = 0 ; irow < nrow ; irow++ ) {
562 val = Zabs(col[2*irow], col[2*irow+1]) ;
563 if ( norm < val ) {
564 norm = val ;
565 }
566 }
567 } else {
568 for ( irow = 0, kk = 0 ; irow < nrow ; irow++, kk += inc1 ) {
569 val = Zabs(col[2*kk], col[2*kk+1]) ;
570 if ( norm < val ) {
571 norm = val ;
572 }
573 }
574 }
575 }
576 return(norm) ; }
577
578 /*--------------------------------------------------------------------*/
579 /*
580 -------------------------------
581 return the one-norm of row irow
582
583 created -- 98apr15, cca
584 -------------------------------
585 */
586 double
A2_oneNormOfRow(A2 * mtx,int irow)587 A2_oneNormOfRow (
588 A2 *mtx,
589 int irow
590 ) {
591 double sum ;
592 double *row ;
593 int inc2, jcol, kk, ncol ;
594 /*
595 ---------------
596 check the input
597 ---------------
598 */
599 if ( mtx == NULL || mtx->entries == NULL
600 || irow < 0 || irow > mtx->n1 ) {
601 fprintf(stderr, "\n fatal error in A2_oneNormOfRow(%p,%d)"
602 "\n bad input\n", mtx, irow) ;
603 exit(-1) ;
604 }
605 if ( ! (A2_IS_REAL(mtx) || A2_IS_COMPLEX(mtx)) ) {
606 fprintf(stderr, "\n fatal error in A2_oneNormOfRow(%p,%d)"
607 "\n bad type %d, must be SPOOLES_REAL or SPOOLES_COMPLEX\n",
608 mtx, irow, mtx->type) ;
609 exit(-1) ;
610 }
611 ncol = mtx->n2 ;
612 sum = 0.0 ;
613 if ( A2_IS_REAL(mtx) ) {
614 row = mtx->entries + irow*mtx->inc1 ;
615 if ( (inc2 = mtx->inc2) == 1 ) {
616 for ( jcol = 0 ; jcol < ncol ; jcol++ ) {
617 sum += fabs(row[jcol]) ;
618 }
619 } else {
620 for ( jcol = 0, kk = 0 ; jcol < ncol ; jcol++, kk += inc2 ) {
621 sum += fabs(row[kk]) ;
622 }
623 }
624 } else if ( A2_IS_COMPLEX(mtx) ) {
625 row = mtx->entries + 2*irow*mtx->inc1 ;
626 if ( (inc2 = mtx->inc2) == 1 ) {
627 for ( jcol = 0 ; jcol < ncol ; jcol++ ) {
628 sum += Zabs(row[2*jcol], row[2*jcol+1]) ;
629 }
630 } else {
631 for ( jcol = 0, kk = 0 ; jcol < ncol ; jcol++, kk += inc2 ) {
632 sum += Zabs(row[2*kk], row[2*kk+1]) ;
633 }
634 }
635 }
636 return(sum) ; }
637
638 /*--------------------------------------------------------------------*/
639 /*
640 -------------------------------
641 return the two-norm of row irow
642
643 created -- 98apr15, cca
644 -------------------------------
645 */
646 double
A2_twoNormOfRow(A2 * mtx,int irow)647 A2_twoNormOfRow (
648 A2 *mtx,
649 int irow
650 ) {
651 double sum ;
652 double *row ;
653 int inc2, jcol, kk, ncol ;
654 /*
655 ---------------
656 check the input
657 ---------------
658 */
659 if ( mtx == NULL || mtx->entries == NULL
660 || irow < 0 || irow > mtx->n1 ) {
661 fprintf(stderr, "\n fatal error in A2_twoNormOfRow(%p,%d)"
662 "\n bad input\n", mtx, irow) ;
663 exit(-1) ;
664 }
665 if ( ! (A2_IS_REAL(mtx) || A2_IS_COMPLEX(mtx)) ) {
666 fprintf(stderr, "\n fatal error in A2_twoNormOfRow(%p,%d)"
667 "\n bad type %d, must be SPOOLES_REAL or SPOOLES_COMPLEX\n",
668 mtx, irow, mtx->type) ;
669 exit(-1) ;
670 }
671 ncol = mtx->n2 ;
672 sum = 0.0 ;
673 if ( A2_IS_REAL(mtx) ) {
674 row = mtx->entries + irow*mtx->inc1 ;
675 if ( (inc2 = mtx->inc2) == 1 ) {
676 for ( jcol = 0 ; jcol < ncol ; jcol++ ) {
677 sum += row[jcol]*row[jcol] ;
678 }
679 } else {
680 for ( jcol = 0, kk = 0 ; jcol < ncol ; jcol++, kk += inc2 ) {
681 sum += row[kk]*row[kk] ;
682 }
683 }
684 } else if ( A2_IS_COMPLEX(mtx) ) {
685 row = mtx->entries + 2*irow*mtx->inc1 ;
686 if ( (inc2 = mtx->inc2) == 1 ) {
687 for ( jcol = 0 ; jcol < ncol ; jcol++ ) {
688 sum += row[2*jcol]*row[2*jcol] + row[2*jcol+1]*row[2*jcol+1] ;
689 }
690 } else {
691 for ( jcol = 0, kk = 0 ; jcol < ncol ; jcol++, kk += inc2 ) {
692 sum += row[2*kk]*row[2*kk] + row[2*kk+1]*row[2*kk+1] ;
693 }
694 }
695 }
696 sum = sqrt(sum) ;
697
698 return(sum) ; }
699
700 /*--------------------------------------------------------------------*/
701 /*
702 ------------------------------------
703 return the infinity-norm of row irow
704
705 created -- 98apr15, cca
706 ------------------------------------
707 */
708 double
A2_infinityNormOfRow(A2 * mtx,int irow)709 A2_infinityNormOfRow (
710 A2 *mtx,
711 int irow
712 ) {
713 double norm, val ;
714 double *row ;
715 int inc2, jcol, kk, ncol ;
716 /*
717 ---------------
718 check the input
719 ---------------
720 */
721 if ( mtx == NULL || mtx->entries == NULL
722 || irow < 0 || irow > mtx->n1 ) {
723 fprintf(stderr, "\n fatal error in A2_infinityNormOfRow(%p,%d)"
724 "\n bad input\n", mtx, irow) ;
725 exit(-1) ;
726 }
727 if ( ! (A2_IS_REAL(mtx) || A2_IS_COMPLEX(mtx)) ) {
728 fprintf(stderr, "\n fatal error in A2_infinityNormOfRow(%p,%d)"
729 "\n bad type %d, must be SPOOLES_REAL or SPOOLES_COMPLEX\n",
730 mtx, irow, mtx->type) ;
731 exit(-1) ;
732 }
733 ncol = mtx->n2 ;
734 norm = 0.0 ;
735 if ( A2_IS_REAL(mtx) ) {
736 row = mtx->entries + irow*mtx->inc1 ;
737 if ( (inc2 = mtx->inc2) == 1 ) {
738 for ( jcol = 0 ; jcol < ncol ; jcol++ ) {
739 val = fabs(row[jcol]) ;
740 if ( norm < val ) {
741 norm = val ;
742 }
743 }
744 } else {
745 for ( jcol = 0, kk = 0 ; jcol < ncol ; jcol++, kk += inc2 ) {
746 val = fabs(row[kk]) ;
747 if ( norm < val ) {
748 norm = val ;
749 }
750 }
751 }
752 } else if ( A2_IS_COMPLEX(mtx) ) {
753 row = mtx->entries + 2*irow*mtx->inc1 ;
754 if ( (inc2 = mtx->inc2) == 1 ) {
755 for ( jcol = 0 ; jcol < ncol ; jcol++ ) {
756 val = Zabs(row[2*jcol], row[2*jcol+1]) ;
757 if ( norm < val ) {
758 norm = val ;
759 }
760 }
761 } else {
762 for ( jcol = 0, kk = 0 ; jcol < ncol ; jcol++, kk += inc2 ) {
763 val = Zabs(row[2*kk], row[2*kk+1]) ;
764 if ( norm < val ) {
765 norm = val ;
766 }
767 }
768 }
769 }
770 return(norm) ; }
771
772 /*--------------------------------------------------------------------*/
773