SpM Handbook 1.2.4
Loading...
Searching...
No Matches
d_spm_norm.c
Go to the documentation of this file.
1/**
2 * @file d_spm_norm.c
3 *
4 * SParse Matrix package norm routine.
5 *
6 * @copyright 2016-2024 Bordeaux INP, CNRS (LaBRI UMR 5800), Inria,
7 * Univ. Bordeaux. All rights reserved.
8 *
9 * @version 1.2.4
10 * @author Mathieu Faverge
11 * @author Pierre Ramet
12 * @author Tony Delarue
13 * @author Matias Hastaran
14 * @author Alycia Lisito
15 * @date 2024-06-25
16 *
17 * @generated from /builds/2mk6rsew/0/fpruvost/spm/src/z_spm_norm.c, normal z -> d, Fri Nov 29 11:34:29 2024
18 *
19 * @ingroup spm_dev_norm
20 * @{
21 *
22 **/
23#include "common.h"
24#include <lapacke.h>
25#include <cblas.h>
26#include "frobeniusupdate.h"
27
28#if !defined(LAPACKE_WITH_LASSQ)
29/**
30 *******************************************************************************
31 *
32 * @brief Updates the values scale and sumsq such that
33 *
34 * ( scale**2 )*sumsq = x( 1 )**2 +...+ x( n )**2 + ( scale**2 )*sumsq,
35 *
36 * This routine is inspired from LAPACK dlassq function.
37 *
38 *******************************************************************************
39 *
40 * @param[in] n
41 * The number of elements in the vector
42 *
43 * @param[in] x
44 * The vector of size abs(n * incx)
45 *
46 * @param[in] incx
47 * The increment between two elements in the vector x.
48 *
49 * @param[inout] scale
50 * On entry, the former scale
51 * On exit, the update scale to take into account the value
52 *
53 * @param[inout] sumsq
54 * On entry, the former sumsq
55 * On exit, the update sumsq to take into account the value
56 *
57 *******************************************************************************
58 *
59 * @return SPM_SUCESS to match the prototype of LAPACKE_dlassq_work
60 *
61 *******************************************************************************/
62static inline int
63__spm_dlassq( spm_int_t n, double *x,
64 spm_int_t incx, double *scale, double *sumsq )
65{
66 spm_int_t i;
67
68 for( i=0; i<n; i++, x+=incx ) {
69#if defined(PRECISION_z) || defined(PRECISION_c)
70 double val;
71 val = ( *x );
72 frobenius_update( 1, scale, sumsq, &val );
73 val = ( *x );
74 frobenius_update( 1, scale, sumsq, &val );
75#else
76 frobenius_update( 1, scale, sumsq, x );
77#endif
78 }
79 return 0;
80}
81
82/**
83 *******************************************************************************
84 *
85 * @brief TODO
86 *
87 *******************************************************************************/
88#define LAPACKE_dlassq_work( _n_, _X_, _incx_, _scale_, _sumsq_ ) \
89 __spm_dlassq( (_n_), (_X_), (_incx_), (_scale_), (_sumsq_) )
90#endif
91
92#if defined(SPM_WITH_MPI)
93/**
94 *******************************************************************************
95 *
96 * @brief MPI reduce operator to merge frobenius partial results together.
97 *
98 *******************************************************************************
99 *
100 * @param[in] dist
101 * TODO
102 *
103 * @param[inout] loc
104 * TODO
105 *
106 * @param[in] len
107 * TODO
108 *
109 * @param[in] dtype
110 * TODO
111 *
112 *******************************************************************************/
113void
114d_spm_frobenius_merge( double *dist,
115 double *loc,
116 int *len,
117 MPI_Datatype *dtype )
118{
119 assert( *len == 2 );
120 frobenius_merge( dist[0], dist[1], loc, loc+1 );
121 (void)len;
122 (void)dtype;
123}
124#endif
125
126/**
127 *******************************************************************************
128 *
129 * @brief Compute the Frobenius norm of a diagonal element within a
130 * symmetric/symmetric matrix with column/row major storage
131 *
132 * Note that column major is using the low triangular part only of the diagonal
133 * element matrices, and row major, by symmetry, is using only the upper
134 * triangular part.
135 *
136 * The comments in the code are made for column major storage.
137 *
138 *******************************************************************************
139 *
140 * @param[in] dofs
141 * TODO
142 *
143 * @param[in] valptr
144 * TODO
145 *
146 * @param[inout] data
147 * TODO
148 *
149 *
150 *******************************************************************************/
151static inline void
152d_spm_frobenius_elt_sym_diag( spm_int_t dofs,
153 const double *valptr,
154 double *data )
155{
156 spm_int_t ii, jj;
157
158 for(jj=0; jj<dofs; jj++)
159 {
160 /* Skip unused upper triangular part */
161 for(ii=0; ii<jj; ii++) {
162 valptr++;
163#if defined(PRECISION_z) || defined(PRECISION_c)
164 valptr++;
165#endif
166 }
167
168 /* Diagonal element */
169 frobenius_update( 1, data, data + 1, valptr );
170#if defined(PRECISION_z) || defined(PRECISION_c)
171 valptr++;
172 frobenius_update( 1, data, data + 1, valptr );
173#endif
174 valptr++;
175
176 for(ii=jj+1; ii<dofs; ii++, valptr++)
177 {
178 frobenius_update( 2, data, data + 1, valptr );
179
180#if defined(PRECISION_z) || defined(PRECISION_c)
181 valptr++;
182 frobenius_update( 2, data, data + 1, valptr );
183#endif
184 }
185 }
186}
187
188/**
189 *******************************************************************************
190 *
191 * @brief Compute the Frobenius norm of an off-diagonal element matrix in the
192 * symmetric/symmetric case
193 *
194 *******************************************************************************
195 *
196 * @param[in] nbelts
197 * TODO
198 *
199 * @param[in] valptr
200 * TODO
201 *
202 * @param[inout] data
203 *
204 *******************************************************************************/
205static inline void
206d_spm_frobenius_elt_sym_offd( spm_int_t nbelts,
207 const double *valptr,
208 double *data )
209{
210 spm_int_t ii;
211
212 for(ii=0; ii<nbelts; ii++, valptr++)
213 {
214 frobenius_update( 2, data, data + 1, valptr );
215
216#if defined(PRECISION_z) || defined(PRECISION_c)
217 valptr++;
218 frobenius_update( 2, data, data + 1, valptr );
219#endif
220 }
221}
222
223/**
224 *******************************************************************************
225 *
226 * @brief Compute the Frobenius norm of any element matrix in the
227 * symmetric/symmetric case
228 *
229 *******************************************************************************
230 *
231 * @param[in] row
232 * TODO
233 *
234 * @param[in] dofi
235 * TODO
236 *
237 * @param[in] col
238 * TODO
239 *
240 * @param[in] dofj
241 * TODO
242 *
243 * @param[in] valptr
244 * TODO
245 *
246 * @param[inout] data
247 * TODO
248 *
249 *
250 *******************************************************************************/
251static inline void
252d_spm_frobenius_elt_sym( spm_int_t row,
253 spm_int_t dofi,
254 spm_int_t col,
255 spm_int_t dofj,
256 const double *valptr,
257 double *data )
258{
259 if ( row == col ) {
260 assert( dofi == dofj );
261 d_spm_frobenius_elt_sym_diag( dofi, (const double*)valptr, data );
262 }
263 else {
264 d_spm_frobenius_elt_sym_offd( dofi * dofj, (const double*)valptr, data );
265 }
266}
267
268/**
269 *******************************************************************************
270 *
271 * @brief Compute the Frobenius norm of a symmetrix/symmetric CSC matrix
272 *
273 *******************************************************************************
274 *
275 * @param[in] spm
276 * The spm from which the norm need to be computed.
277 *
278 * @param[in,out] data
279 * TODO
280 *
281 *******************************************************************************/
282static inline void
283d_spmFrobeniusNorm_csc( const spmatrix_t *spm,
284 double *data )
285{
286 spm_int_t j, k, baseval;
287 spm_int_t ig, dofi, row;
288 spm_int_t jg, dofj, col;
289 const spm_int_t *colptr;
290 const spm_int_t *rowptr;
291 const spm_int_t *dofs;
292 const spm_int_t *loc2glob;
293 const double *valptr;
294
295 assert( spm->fmttype == SpmCSC );
296 assert( spm->flttype == SpmDouble );
297
298 baseval = spm->baseval;
299 colptr = spm->colptr;
300 rowptr = spm->rowptr;
301 valptr = (double*)(spm->values);
302 dofs = spm->dofs;
303 loc2glob = spm->loc2glob;
304
305 for(j=0; j<spm->n; j++, colptr++, loc2glob++)
306 {
307 jg = spm->replicated ? j : (*loc2glob) - baseval;
308 if ( spm->dof > 0 ) {
309 dofj = spm->dof;
310 col = spm->dof * jg;
311 }
312 else {
313 dofj = dofs[jg+1] - dofs[jg];
314 col = dofs[jg] - baseval;
315 }
316
317 for(k=colptr[0]; k<colptr[1]; k++, rowptr++)
318 {
319 ig = (*rowptr - baseval);
320 if ( spm->dof > 0 ) {
321 dofi = spm->dof;
322 row = spm->dof * ig;
323 }
324 else {
325 dofi = dofs[ig+1] - dofs[ig];
326 row = dofs[ig] - baseval;
327 }
328
329 d_spm_frobenius_elt_sym( row, dofi, col, dofj, valptr, data );
330 valptr += dofi * dofj;
331 }
332 }
333}
334
335/**
336 *******************************************************************************
337 *
338 * @brief Compute the Frobenius norm of a symmetrix/symmetric CSR matrix
339 *
340 *******************************************************************************
341 *
342 * @param[in] spm
343 * The spm from which the norm need to be computed.
344 *
345 * @param[in,out] data
346 * TODO
347 *
348 *******************************************************************************/
349static inline void
350d_spmFrobeniusNorm_csr( const spmatrix_t *spm,
351 double *data )
352{
353 spm_int_t i, k, baseval;
354 spm_int_t ig, dofi, row;
355 spm_int_t jg, dofj, col;
356 const spm_int_t *colptr;
357 const spm_int_t *rowptr;
358 const spm_int_t *dofs;
359 const spm_int_t *loc2glob;
360 const double *valptr;
361
362 assert( spm->fmttype == SpmCSR );
363 assert( spm->flttype == SpmDouble );
364
365 baseval = spm->baseval;
366
367 colptr = spm->colptr;
368 rowptr = spm->rowptr;
369 valptr = (double*)(spm->values);
370 dofs = spm->dofs;
371 loc2glob = spm->loc2glob;
372
373 for(i=0; i<spm->n; i++, rowptr++, loc2glob++)
374 {
375 ig = spm->replicated ? i : (*loc2glob) - baseval;
376 if ( spm->dof > 0 ) {
377 dofi = spm->dof;
378 row = spm->dof * ig;
379 }
380 else {
381 dofi = dofs[ig+1] - dofs[ig];
382 row = dofs[ig] - baseval;
383 }
384
385 for(k=rowptr[0]; k<rowptr[1]; k++, colptr++)
386 {
387 jg = (*colptr - baseval);
388 if ( spm->dof > 0 ) {
389 dofj = spm->dof;
390 col = spm->dof * jg;
391 }
392 else {
393 dofj = dofs[jg+1] - dofs[jg];
394 col = dofs[jg] - baseval;
395 }
396
397 d_spm_frobenius_elt_sym( row, dofi, col, dofj, valptr, data );
398 valptr += dofi * dofj;
399 }
400 }
401}
402
403/**
404 *******************************************************************************
405 *
406 * @brief Compute the Frobenius norm of a symmetrix/symmetric IJV matrix
407 *
408 *******************************************************************************
409 *
410 * @param[in] spm
411 * The spm from which the norm need to be computed.
412 *
413 * @param[in,out] data
414 * TODO
415 *
416 *******************************************************************************/
417static inline void
418d_spmFrobeniusNorm_ijv( const spmatrix_t *spm,
419 double *data )
420{
421 spm_int_t k, baseval;
422 spm_int_t i, dofi, row;
423 spm_int_t j, dofj, col;
424 const spm_int_t *colptr;
425 const spm_int_t *rowptr;
426 const spm_int_t *dofs;
427 const double *valptr;
428
429 assert( spm->fmttype == SpmIJV );
430 assert( spm->flttype == SpmDouble );
431
432 baseval = spm->baseval;
433
434 colptr = spm->colptr;
435 rowptr = spm->rowptr;
436 valptr = (double*)(spm->values);
437 dofs = spm->dofs;
438
439 for(k=0; k<spm->nnz; k++, rowptr++, colptr++)
440 {
441 i = *rowptr - baseval;
442 j = *colptr - baseval;
443
444 if ( spm->dof > 0 ) {
445 dofi = spm->dof;
446 row = spm->dof * i;
447 dofj = spm->dof;
448 col = spm->dof * j;
449 }
450 else {
451 dofi = dofs[i+1] - dofs[i];
452 row = dofs[i] - baseval;
453 dofj = dofs[j+1] - dofs[j];
454 col = dofs[j] - baseval;
455 }
456
457 d_spm_frobenius_elt_sym( row, dofi, col, dofj, valptr, data );
458 valptr += dofi * dofj;
459 }
460}
461
462/**
463 *******************************************************************************
464 *
465 * @brief Compute the Frobenius norm of the given spm structure.
466 *
467 * ||A|| = sqrt( sum( a_ij ^ 2 ) )
468 *
469 *******************************************************************************
470 *
471 * @param[in] spm
472 * The spm from which the norm need to be computed.
473 *
474 *******************************************************************************
475 *
476 * @return The computed frobenius norm
477 *
478 *******************************************************************************/
479double
480d_spmFrobeniusNorm( const spmatrix_t *spm )
481{
482 double data[] = { 0., 1. }; /* Scale, Sum */
483
484 if (spm->mtxtype == SpmGeneral) {
485 const double *valptr = (double*)spm->values;
486 spm_int_t i;
487
488 for(i=0; i <spm->nnzexp; i++, valptr++) {
489 frobenius_update( 1, data, data + 1, valptr );
490
491#if defined(PRECISION_z) || defined(PRECISION_c)
492 valptr++;
493 frobenius_update( 1, data, data + 1, valptr );
494#endif
495 }
496 }
497 else {
498 switch( spm->fmttype ) {
499 case SpmCSC:
500 d_spmFrobeniusNorm_csc( spm, data );
501 break;
502
503 case SpmCSR:
504 d_spmFrobeniusNorm_csr( spm, data );
505 break;
506
507 case SpmIJV:
508 default:
509 d_spmFrobeniusNorm_ijv( spm, data );
510 }
511 }
512
513#if defined(SPM_WITH_MPI)
514 if ( !(spm->replicated) && (spm->clustnbr > 1) ) {
515 MPI_Op merge;
516 MPI_Op_create( (MPI_User_function *)d_spm_frobenius_merge, 1, &merge );
517 MPI_Allreduce( MPI_IN_PLACE, data, 2, SPM_MPI_DOUBLE, merge, spm->comm );
518 MPI_Op_free( &merge );
519 }
520#endif
521
522 return data[0] * sqrt( data[1] );
523}
524
525/**
526 *******************************************************************************
527 *
528 * @brief Compute the Max norm of the given spm structure.
529 *
530 * ||A|| = max( abs(a_ij) )
531 *
532 *******************************************************************************
533 *
534 * @param[in] spm
535 * The spm from which the norm need to be computed.
536 *
537 *******************************************************************************
538 *
539 * @return The computed max norm
540 *
541 *******************************************************************************/
542double
543d_spmMaxNorm( const spmatrix_t *spm )
544{
545 spm_int_t i;
546 const double *valptr = (double *)spm->values;
547 double tmp, norm = 0.;
548
549 for(i=0; i <spm->nnzexp; i++, valptr++) {
550 tmp = fabs( *valptr );
551 norm = (norm > tmp) ? norm : tmp;
552 }
553
554#if defined(SPM_WITH_MPI)
555 if ( !(spm->replicated) && (spm->clustnbr > 1) ) {
556 MPI_Allreduce( MPI_IN_PLACE, &norm, 1, MPI_DOUBLE, MPI_MAX, spm->comm );
557 }
558#endif
559
560 return norm;
561}
562
563/**
564 *******************************************************************************
565 *
566 * @brief Compute the sum array for the one/inf norms of a diagonal element
567 * within a symmetric/symmetric matrix with column/row major storage.
568 *
569 * Note that column major is using the low triangular part only of the diagonal
570 * element matrices, and row major, by symmetry, is using only the upper
571 * triangular part.
572 *
573 * The comments in the code are made for column major storage.
574 *
575 *******************************************************************************
576 *
577 * @param[in] row
578 * TODO
579 *
580 * @param[in] dofi
581 * TODO
582 *
583 * @param[in] valptr
584 * TODO
585 *
586 * @param[inout] sumtab
587 * TODO
588 *
589 *******************************************************************************/
590static inline void
591d_spm_oneinf_elt_sym_diag( spm_int_t row,
592 spm_int_t dofi,
593 const double *valptr,
594 double *sumtab )
595{
596 spm_int_t ii, jj;
597
598 sumtab += row;
599
600 for(jj=0; jj<dofi; jj++)
601 {
602 /* Skip unused upper triangular part */
603 for(ii=0; ii<jj; ii++) {
604 valptr++;
605 }
606
607 /* Diagonal element */
608 sumtab[jj] += fabs( *valptr );
609 valptr++;
610
611 for(ii=jj+1; ii<dofi; ii++, valptr++)
612 {
613 /* Lower part */
614 sumtab[ii] += fabs( *valptr );
615 /* Upper part */
616 sumtab[jj] += fabs( *valptr );
617 }
618 }
619}
620
621/**
622 *******************************************************************************
623 *
624 * @brief Compute the sum array for the one/inf norms of a general element.
625 *
626 * We can observe two cases A and B;
627 * ________ _ ________ _
628 * | | | | | | |________| | |
629 * | | | | |A| OR |________| |B|
630 * |__|__|__| |_| |________| |_|
631 * ________ ________
632 * |___B____| |___A____|
633 *
634 * | One Norm | Inf norm |
635 * -------------+----------+----------+
636 * Column Major | B | A |
637 * -------------+----------+----------+
638 * Row Major | A | B |
639 * -------------+----------+----------+
640 *
641 * @warning: The sumtab must be shifted at the right place on input
642 *
643 *******************************************************************************
644 *
645 * @param[in] dofi
646 * TODO
647 *
648 * @param[in] dofj
649 * TODO
650 *
651 * @param[in] valptr
652 * TODO
653 *
654 * @param[inout] sumtab
655 * TODO
656 *
657 *******************************************************************************/
658static inline void
659d_spm_oneinf_elt_gen_A( spm_int_t dofi,
660 spm_int_t dofj,
661 const double *valptr,
662 double *sumtab )
663{
664 spm_int_t ii, jj;
665
666 for(jj=0; jj<dofj; jj++)
667 {
668 for(ii=0; ii<dofi; ii++, valptr++)
669 {
670 sumtab[ii] += fabs( *valptr );
671 }
672 }
673}
674
675/**
676 *******************************************************************************
677 *
678 * @brief Compute the sum array for the one/inf norms of a general element.
679 *
680 * See d_spm_oneinf_elt_gen_A()
681 *
682 *******************************************************************************
683 *
684 * @param[in] dofi
685 * TODO
686 *
687 * @param[in] dofj
688 * TODO
689 *
690 * @param[in] valptr
691 * TODO
692 *
693 * @param[inout] sumtab
694 * TODO
695 *
696 *******************************************************************************/
697static inline void
698d_spm_oneinf_elt_gen_B( spm_int_t dofi,
699 spm_int_t dofj,
700 const double *valptr,
701 double *sumtab )
702{
703 spm_int_t ii, jj;
704
705 for(jj=0; jj<dofj; jj++, sumtab++)
706 {
707 for(ii=0; ii<dofi; ii++, valptr++)
708 {
709 *sumtab += fabs( *valptr );
710 }
711 }
712}
713
714/**
715 *******************************************************************************
716 *
717 * @brief Compute the sum array for both the one and inf norms for the
718 * off-diagonal elements of symmetric/symmetric element matrices.
719 *
720 * See d_spm_oneinf_elt_gen_A()
721 *
722 *******************************************************************************
723 *
724 * @param[in] row
725 * TODO
726 *
727 * @param[in] dofi
728 * TODO
729 *
730 * @param[in] col
731 * TODO
732 *
733 * @param[in] dofj
734 * TODO
735 *
736 * @param[in] valptr
737 * TODO
738 *
739 * @param[inout] sumtab
740 * TODO
741 *
742 *******************************************************************************/
743static inline void
744d_spm_oneinf_elt_gen_AB( spm_int_t row,
745 spm_int_t dofi,
746 spm_int_t col,
747 spm_int_t dofj,
748 const double *valptr,
749 double *sumtab )
750{
751 double *sumrow = sumtab + row;
752 double *sumcol = sumtab + col;
753 spm_int_t ii, jj;
754
755 for(jj=0; jj<dofj; jj++, sumcol++)
756 {
757 for(ii=0; ii<dofi; ii++, valptr++)
758 {
759 double v = fabs( *valptr );
760 sumrow[ii] += v;
761 *sumcol += v;
762 }
763 }
764}
765
766/**
767 *******************************************************************************
768 *
769 * @brief Compute the sum array for the one and inf norms of a general element.
770 *
771 *******************************************************************************
772 *
773 * @param[in] layout
774 * TODO
775 *
776 * @param[in] row
777 * TODO
778 *
779 * @param[in] dofi
780 * TODO
781 *
782 * @param[in] col
783 * TODO
784 *
785 * @param[in] dofj
786 * TODO
787 *
788 * @param[in] valptr
789 * TODO
790 *
791 * @param[in] ntype
792 * TODO
793 *
794 * @param[inout] sumtab
795 * TODO
796 *
797 *******************************************************************************/
798static inline void
799d_spm_oneinf_elt_gen( spm_layout_t layout,
800 spm_int_t row,
801 spm_int_t dofi,
802 spm_int_t col,
803 spm_int_t dofj,
804 const double *valptr,
805 spm_normtype_t ntype,
806 double *sumtab )
807{
808 if ( layout == SpmColMajor ) {
809 if ( ntype == SpmInfNorm ) {
810 d_spm_oneinf_elt_gen_A( dofi, dofj, valptr, sumtab + row );
811 }
812 else {
813 assert( ntype == SpmOneNorm );
814 d_spm_oneinf_elt_gen_B( dofi, dofj, valptr, sumtab + col );
815 }
816 }
817 else {
818 if ( ntype == SpmInfNorm ) {
819 d_spm_oneinf_elt_gen_B( dofj, dofi, valptr, sumtab + row );
820 }
821 else {
822 assert( ntype == SpmOneNorm );
823 d_spm_oneinf_elt_gen_A( dofj, dofi, valptr, sumtab + col );
824 }
825 }
826}
827
828/**
829 *******************************************************************************
830 *
831 * @brief Compute the sum array for both the one and inf norms for the
832 * off-diagonal elements of symmetric/symmetric element matrices in either
833 * column or row major layout.
834 *
835 *******************************************************************************
836 *
837 * @param[in] layout
838 * TODO
839 *
840 * @param[in] row
841 * TODO
842 *
843 * @param[in] dofi
844 * TODO
845 *
846 * @param[in] col
847 * TODO
848 *
849 * @param[in] dofj
850 * TODO
851 *
852 * @param[in] valptr
853 * TODO
854 *
855 * @param[inout] sumtab
856 * TODO
857 *
858 *******************************************************************************/
859static inline void
860d_spm_oneinf_elt_sym_offd( spm_layout_t layout,
861 spm_int_t row,
862 spm_int_t dofi,
863 spm_int_t col,
864 spm_int_t dofj,
865 const double *valptr,
866 double *sumtab )
867{
868 if ( layout == SpmColMajor ) {
869 d_spm_oneinf_elt_gen_AB( row, dofi, col, dofj, valptr, sumtab );
870 }
871 else {
872 d_spm_oneinf_elt_gen_AB( col, dofj, row, dofi, valptr, sumtab );
873 }
874}
875
876/**
877 *******************************************************************************
878 *
879 * @brief Compute the sum array for the one/inf norm for an element matrix.
880 *
881 *******************************************************************************
882 *
883 * @param[in] mtxtype
884 * TODO
885 *
886 * @param[in] layout
887 * TODO
888 *
889 * @param[in] row
890 * TODO
891 *
892 * @param[in] dofi
893 * TODO
894 *
895 * @param[in] col
896 * TODO
897 *
898 * @param[in] dofj
899 * TODO
900 *
901 * @param[in] valptr
902 * TODO
903 *
904 * @param[in] ntype
905 * TODO
906 *
907 * @param[inout] sumtab
908 * TODO
909 *
910 *******************************************************************************/
911static inline void
912d_spm_oneinf_elt( spm_mtxtype_t mtxtype,
913 spm_layout_t layout,
914 spm_int_t row,
915 spm_int_t dofi,
916 spm_int_t col,
917 spm_int_t dofj,
918 const double *valptr,
919 spm_normtype_t ntype,
920 double *sumtab )
921{
922 if ( mtxtype == SpmGeneral ) {
923 d_spm_oneinf_elt_gen( layout, row, dofi, col, dofj, valptr, ntype, sumtab );
924 }
925 else {
926 if ( row == col ) {
927 d_spm_oneinf_elt_sym_diag( row, dofi, valptr, sumtab );
928 }
929 else {
930 d_spm_oneinf_elt_sym_offd( layout, row, dofi, col, dofj, valptr, sumtab );
931 }
932 }
933}
934
935/**
936 *******************************************************************************
937 *
938 * @brief Compute the one/inf norm of an spm CSC structure.
939 *
940 *******************************************************************************
941 *
942 * @param[in] ntype
943 * TODO
944 *
945 * @param[in] spm
946 * TODO
947 *
948 * @param[inout] sumtab
949 * TODO
950 *
951 *******************************************************************************/
952static inline void
953d_spmOneInfNorm_csc( spm_normtype_t ntype,
954 const spmatrix_t *spm,
955 double *sumtab )
956{
957 spm_int_t i, j, ig, jg, col, row;
958 spm_int_t dofi, dofj, dof, baseval;
959 const spm_int_t *colptr, *rowptr, *loc2glob, *dofs;
960 const double *valptr;
961
962 baseval = spm->baseval;
963 colptr = spm->colptr;
964 rowptr = spm->rowptr;
965 valptr = (const double *)(spm->values);
966 loc2glob = spm->loc2glob;
967 dofs = spm->dofs;
968 dof = spm->dof;
969 for(j=0; j<spm->n; j++, colptr++, loc2glob++)
970 {
971 jg = spm->replicated ? j : (*loc2glob) - baseval;
972 if ( dof > 0 ) {
973 dofj = dof;
974 col = dof * jg;
975 }
976 else {
977 dofj = dofs[jg+1] - dofs[jg];
978 col = dofs[jg] - baseval;
979 }
980
981 for(i=colptr[0]; i<colptr[1]; i++, rowptr++)
982 {
983 ig = (*rowptr - baseval);
984 if ( dof > 0 ) {
985 dofi = dof;
986 row = dof * ig;
987 }
988 else {
989 dofi = dofs[ig+1] - dofs[ig];
990 row = dofs[ig] - baseval;
991 }
992
993 d_spm_oneinf_elt( spm->mtxtype, spm->layout,
994 row, dofi, col, dofj, valptr,
995 ntype, sumtab );
996 valptr += dofi * dofj;
997 }
998 }
999}
1000
1001/**
1002 *******************************************************************************
1003 *
1004 * @brief Compute the one/inf norm of an spm CSR structure.
1005 *
1006 *******************************************************************************
1007 *
1008 * @param[in] ntype
1009 * TODO
1010 *
1011 * @param[in] spm
1012 * TODO
1013 *
1014 * @param[inout] sumtab
1015 * TODO
1016 *
1017 *******************************************************************************/
1018static inline void
1019d_spmOneInfNorm_csr( spm_normtype_t ntype,
1020 const spmatrix_t *spm,
1021 double *sumtab )
1022{
1023 spm_int_t i, j, ig, jg, col, row;
1024 spm_int_t dofi, dofj, dof, baseval;
1025 const spm_int_t *colptr, *rowptr, *loc2glob, *dofs;
1026 const double *valptr;
1027
1028 baseval = spm->baseval;
1029 colptr = spm->colptr;
1030 rowptr = spm->rowptr;
1031 valptr = (const double *)(spm->values);
1032 loc2glob = spm->loc2glob;
1033 dofs = spm->dofs;
1034 dof = spm->dof;
1035 for(i=0; i<spm->n; i++, rowptr++, loc2glob++)
1036 {
1037 ig = spm->replicated ? i : (*loc2glob) - baseval;
1038 if ( dof > 0 ) {
1039 dofi = dof;
1040 row = dof * ig;
1041 }
1042 else {
1043 dofi = dofs[ig+1] - dofs[ig];
1044 row = dofs[ig] - baseval;
1045 }
1046
1047 for(j=rowptr[0]; j<rowptr[1]; j++, colptr++)
1048 {
1049 jg = (*colptr - baseval);
1050 if ( dof > 0 ) {
1051 dofj = dof;
1052 col = dof * jg;
1053 }
1054 else {
1055 dofj = dofs[jg+1] - dofs[jg];
1056 col = dofs[jg] - baseval;
1057 }
1058
1059 d_spm_oneinf_elt( spm->mtxtype, spm->layout,
1060 row, dofi, col, dofj, valptr,
1061 ntype, sumtab );
1062 valptr += dofi * dofj;
1063 }
1064 }
1065}
1066
1067/**
1068 *******************************************************************************
1069 *
1070 * @brief Compute the one/inf norm of an spm IJV structure.
1071 *
1072 *******************************************************************************
1073 *
1074 * @param[in] ntype
1075 * TODO
1076 *
1077 * @param[in] spm
1078 * TODO
1079 *
1080 * @param[inout] sumtab
1081 * TODO
1082 *
1083 *******************************************************************************/
1084static inline void
1085d_spmOneInfNorm_ijv( spm_normtype_t ntype,
1086 const spmatrix_t *spm,
1087 double *sumtab )
1088{
1089 spm_int_t k, ig, jg, col, row;
1090 spm_int_t dofi, dofj, dof, baseval;
1091 spm_int_t *colptr, *rowptr, *dofs;
1092 double *valptr;
1093
1094 baseval = spm->baseval;
1095 colptr = spm->colptr;
1096 rowptr = spm->rowptr;
1097 valptr = (double *)(spm->values);
1098 dofs = spm->dofs;
1099 dof = spm->dof;
1100
1101 for(k=0; k<spm->nnz; k++, rowptr++, colptr++)
1102 {
1103 ig = *rowptr - baseval;
1104 jg = *colptr - baseval;
1105
1106 if ( dof > 0 ) {
1107 dofi = dof;
1108 row = dof * ig;
1109 dofj = dof;
1110 col = dof * jg;
1111 }
1112 else {
1113 dofi = dofs[ig+1] - dofs[ig];
1114 row = dofs[ig] - baseval;
1115 dofj = dofs[jg+1] - dofs[jg];
1116 col = dofs[jg] - baseval;
1117 }
1118
1119 d_spm_oneinf_elt( spm->mtxtype, spm->layout,
1120 row, dofi, col, dofj, valptr,
1121 ntype, sumtab );
1122 valptr += dofi * dofj;
1123 }
1124}
1125
1126/**
1127 *******************************************************************************
1128 *
1129 * @brief Compute the one/inf norm of the given spm structure given by
1130 * the maximum row sum
1131 *
1132 * * SpmOneNorm: ||A|| = max_j( sum_i(|a_ij|) )
1133 * * SpmInfNorm: ||A|| = max_i( sum_j(|a_ij|) )
1134 *
1135 *******************************************************************************
1136 *
1137 * @param[in] ntype
1138 * The type of norm to compute.
1139 *
1140 * @param[in] spm
1141 * The spm from which the norm need to be computed.
1142 *
1143 *******************************************************************************
1144 *
1145 * @return The computed one norm
1146 *
1147 *******************************************************************************/
1148static inline double
1149d_spmOneInfNorm( spm_normtype_t ntype,
1150 const spmatrix_t *spm )
1151{
1152 spm_int_t k;
1153 double *sumtab = calloc( spm->gNexp, sizeof(double) );
1154 double norm = 0.;
1155
1156 switch( spm->fmttype ) {
1157 case SpmCSC:
1158 d_spmOneInfNorm_csc( ntype, spm, sumtab );
1159 break;
1160
1161 case SpmCSR:
1162 d_spmOneInfNorm_csr( ntype, spm, sumtab );
1163 break;
1164
1165 case SpmIJV:
1166 default:
1167 d_spmOneInfNorm_ijv( ntype, spm, sumtab );
1168 }
1169
1170#if defined(SPM_WITH_MPI)
1171 if ( !(spm->replicated) && (spm->clustnbr > 1) ) {
1172 MPI_Allreduce( MPI_IN_PLACE, sumtab, spm->gNexp, MPI_DOUBLE, MPI_SUM, spm->comm );
1173 }
1174#endif
1175
1176 /* Look for the maximum */
1177 {
1178 const double *sumtmp = sumtab;
1179 for( k=0; k<spm->gNexp; k++, sumtmp++ )
1180 {
1181 if( norm < *sumtmp ) {
1182 norm = *sumtmp;
1183 }
1184 }
1185 }
1186
1187 free( sumtab );
1188 return norm;
1189}
1190
1191/**
1192 *******************************************************************************
1193 *
1194 * @brief Compute the norm of an spm matrix
1195 *
1196 *******************************************************************************
1197 *
1198 * @param[in] ntype
1199 * = SpmMaxNorm: Max norm
1200 * = SpmOneNorm: One norm
1201 * = SpmInfNorm: Infinity norm
1202 * = SpmFrobeniusNorm: Frobenius norm
1203 *
1204 * @param[in] spm
1205 * The spm structure describing the matrix.
1206 *
1207 *******************************************************************************
1208 *
1209 * @return The norm of the spm matrix
1210 * -1 when error occurs or with pattern only
1211 *
1212 *******************************************************************************/
1213double
1214d_spmNorm( spm_normtype_t ntype,
1215 const spmatrix_t *spm )
1216{
1217 double norm = 0.;
1218
1219 if ( spm == NULL ) {
1220 return -1.;
1221 }
1222
1223 switch( ntype ) {
1224 case SpmMaxNorm:
1225 norm = d_spmMaxNorm( spm );
1226 break;
1227
1228 case SpmInfNorm:
1229 case SpmOneNorm:
1230 norm = d_spmOneInfNorm( ntype, spm );
1231 break;
1232
1233 case SpmFrobeniusNorm:
1234 norm = d_spmFrobeniusNorm( spm );
1235 break;
1236
1237 default:
1238 fprintf(stderr, "d_spmNorm: invalid norm type\n");
1239 return -1.;
1240 }
1241
1242 return norm;
1243}
1244
1245/**
1246 *******************************************************************************
1247 *
1248 * @brief Compute the norm of a dense matrix that follows the distribution of an
1249 * spm matrix
1250 *
1251 *******************************************************************************
1252 *
1253 * @param[in] ntype
1254 * = SpmMaxNorm: Max norm
1255 * = SpmOneNorm: One norm
1256 * = SpmInfNorm: Infinity norm
1257 * = SpmFrobeniusNorm: Frobenius norm
1258 *
1259 * @param[in] spm
1260 * The spm structure describing the matrix.
1261 *
1262 * @param[in] n
1263 * The number of columns of the matrix A.
1264 *
1265 * @param[in] A
1266 * The matrix A of size lda-by-n.
1267 *
1268 * @param[in] lda
1269 * The leading dimension of the matrix A. Must be >= max(1, spm->nexp).
1270 *
1271 *******************************************************************************
1272 *
1273 * @return The norm of the spm matrix
1274 * -1 when error occurs or with pattern only
1275 *
1276 *******************************************************************************/
1277double
1278d_spmNormMat( spm_normtype_t ntype,
1279 const spmatrix_t *spm,
1280 spm_int_t n,
1281 const double *A,
1282 spm_int_t lda )
1283{
1284 double norm = 0.;
1285 int j;
1286
1287 if ( spm == NULL ) {
1288 return -1.;
1289 }
1290
1291 switch( ntype ) {
1292 case SpmMaxNorm:
1293 case SpmInfNorm:
1294 norm = LAPACKE_dlange( LAPACK_COL_MAJOR,
1295 ntype == SpmMaxNorm ? 'M' : 'I',
1296 spm->nexp, n, A, lda );
1297#if defined(SPM_WITH_MPI)
1298 if ( !(spm->replicated) && (spm->clustnbr > 1) ) {
1299 MPI_Allreduce( MPI_IN_PLACE, &norm, 1, MPI_DOUBLE,
1300 MPI_MAX, spm->comm );
1301 }
1302#endif
1303 break;
1304
1305 case SpmOneNorm:
1306 {
1307 double *sumtmp;
1308 double *sumtab = calloc( n, sizeof(double) );
1309
1310 sumtmp = sumtab;
1311 for( j=0; j<n; j++, sumtmp++ )
1312 {
1313 *sumtmp = cblas_dasum( spm->nexp, A + j * lda, 1 );
1314 }
1315
1316#if defined(SPM_WITH_MPI)
1317 if ( !(spm->replicated) && (spm->clustnbr > 1) ) {
1318 MPI_Allreduce( MPI_IN_PLACE, sumtab, n, MPI_DOUBLE,
1319 MPI_SUM, spm->comm );
1320 }
1321#endif
1322
1323 /* Look for the maximum */
1324 sumtmp = sumtab;
1325 for( j=0; j<n; j++, sumtmp++ )
1326 {
1327 if( norm < *sumtmp ) {
1328 norm = *sumtmp;
1329 }
1330 }
1331 free( sumtab );
1332 }
1333 break;
1334
1335 case SpmFrobeniusNorm:
1336 {
1337 double data[] = { 0., 1. }; /* Scale, Sum */
1338
1339 for ( j=0; j<n; j++ ) {
1340 /* LAPACKE interface is incorrect and do not have the const yet */
1341 LAPACKE_dlassq_work( spm->nexp, (double*)(A + j * lda), 1, data, data + 1 );
1342 }
1343
1344#if defined(SPM_WITH_MPI)
1345 if ( !(spm->replicated) && (spm->clustnbr > 1) ) {
1346 MPI_Op merge;
1347 MPI_Op_create( (MPI_User_function *)d_spm_frobenius_merge, 1, &merge );
1348 MPI_Allreduce( MPI_IN_PLACE, data, 2, SPM_MPI_DOUBLE, merge, spm->comm );
1349 MPI_Op_free( &merge );
1350 }
1351#endif
1352
1353 norm = data[0] * sqrt( data[1] );
1354 }
1355 break;
1356
1357 default:
1358 fprintf(stderr, "d_spmNorm: invalid norm type\n");
1359 return -1.;
1360 }
1361
1362 return norm;
1363}
1364/**
1365 * @}
1366 */
d_spmOneInfNorm_csc
static void d_spmOneInfNorm_csc(spm_normtype_t ntype, const spmatrix_t *spm, double *sumtab)
Compute the one/inf norm of an spm CSC structure.
Definition d_spm_norm.c:953
d_spm_oneinf_elt_gen_AB
static void d_spm_oneinf_elt_gen_AB(spm_int_t row, spm_int_t dofi, spm_int_t col, spm_int_t dofj, const double *valptr, double *sumtab)
Compute the sum array for both the one and inf norms for the off-diagonal elements of symmetric/symme...
Definition d_spm_norm.c:744
d_spm_oneinf_elt_gen_B
static void d_spm_oneinf_elt_gen_B(spm_int_t dofi, spm_int_t dofj, const double *valptr, double *sumtab)
Compute the sum array for the one/inf norms of a general element.
Definition d_spm_norm.c:698
d_spmMaxNorm
double d_spmMaxNorm(const spmatrix_t *spm)
Compute the Max norm of the given spm structure.
Definition d_spm_norm.c:543
d_spmOneInfNorm_ijv
static void d_spmOneInfNorm_ijv(spm_normtype_t ntype, const spmatrix_t *spm, double *sumtab)
Compute the one/inf norm of an spm IJV structure.
Definition d_spm_norm.c:1085
d_spm_oneinf_elt_sym_diag
static void d_spm_oneinf_elt_sym_diag(spm_int_t row, spm_int_t dofi, const double *valptr, double *sumtab)
Compute the sum array for the one/inf norms of a diagonal element within a symmetric/symmetric matrix...
Definition d_spm_norm.c:591
LAPACKE_dlassq_work
#define LAPACKE_dlassq_work(_n_, _X_, _incx_, _scale_, _sumsq_)
TODO.
Definition d_spm_norm.c:88
d_spmFrobeniusNorm_csc
static void d_spmFrobeniusNorm_csc(const spmatrix_t *spm, double *data)
Compute the Frobenius norm of a symmetrix/symmetric CSC matrix.
Definition d_spm_norm.c:283
d_spmOneInfNorm
static double d_spmOneInfNorm(spm_normtype_t ntype, const spmatrix_t *spm)
Compute the one/inf norm of the given spm structure given by the maximum row sum.
Definition d_spm_norm.c:1149
d_spmFrobeniusNorm
double d_spmFrobeniusNorm(const spmatrix_t *spm)
Compute the Frobenius norm of the given spm structure.
Definition d_spm_norm.c:480
d_spm_oneinf_elt_sym_offd
static void d_spm_oneinf_elt_sym_offd(spm_layout_t layout, spm_int_t row, spm_int_t dofi, spm_int_t col, spm_int_t dofj, const double *valptr, double *sumtab)
Compute the sum array for both the one and inf norms for the off-diagonal elements of symmetric/symme...
Definition d_spm_norm.c:860
d_spm_oneinf_elt_gen_A
static void d_spm_oneinf_elt_gen_A(spm_int_t dofi, spm_int_t dofj, const double *valptr, double *sumtab)
Compute the sum array for the one/inf norms of a general element.
Definition d_spm_norm.c:659
__spm_dlassq
static int __spm_dlassq(spm_int_t n, double *x, spm_int_t incx, double *scale, double *sumsq)
Updates the values scale and sumsq such that.
Definition d_spm_norm.c:63
d_spm_frobenius_elt_sym_diag
static void d_spm_frobenius_elt_sym_diag(spm_int_t dofs, const double *valptr, double *data)
Compute the Frobenius norm of a diagonal element within a symmetric/symmetric matrix with column/row ...
Definition d_spm_norm.c:152
d_spmOneInfNorm_csr
static void d_spmOneInfNorm_csr(spm_normtype_t ntype, const spmatrix_t *spm, double *sumtab)
Compute the one/inf norm of an spm CSR structure.
Definition d_spm_norm.c:1019
d_spmFrobeniusNorm_ijv
static void d_spmFrobeniusNorm_ijv(const spmatrix_t *spm, double *data)
Compute the Frobenius norm of a symmetrix/symmetric IJV matrix.
Definition d_spm_norm.c:418
d_spm_frobenius_elt_sym_offd
static void d_spm_frobenius_elt_sym_offd(spm_int_t nbelts, const double *valptr, double *data)
Compute the Frobenius norm of an off-diagonal element matrix in the symmetric/symmetric case.
Definition d_spm_norm.c:206
d_spm_frobenius_elt_sym
static void d_spm_frobenius_elt_sym(spm_int_t row, spm_int_t dofi, spm_int_t col, spm_int_t dofj, const double *valptr, double *data)
Compute the Frobenius norm of any element matrix in the symmetric/symmetric case.
Definition d_spm_norm.c:252
d_spm_oneinf_elt
static void d_spm_oneinf_elt(spm_mtxtype_t mtxtype, spm_layout_t layout, spm_int_t row, spm_int_t dofi, spm_int_t col, spm_int_t dofj, const double *valptr, spm_normtype_t ntype, double *sumtab)
Compute the sum array for the one/inf norm for an element matrix.
Definition d_spm_norm.c:912
d_spmFrobeniusNorm_csr
static void d_spmFrobeniusNorm_csr(const spmatrix_t *spm, double *data)
Compute the Frobenius norm of a symmetrix/symmetric CSR matrix.
Definition d_spm_norm.c:350
d_spm_oneinf_elt_gen
static void d_spm_oneinf_elt_gen(spm_layout_t layout, spm_int_t row, spm_int_t dofi, spm_int_t col, spm_int_t dofj, const double *valptr, spm_normtype_t ntype, double *sumtab)
Compute the sum array for the one and inf norms of a general element.
Definition d_spm_norm.c:799
spm_layout_t
enum spm_layout_e spm_layout_t
Direction of the matrix storage.
spm_normtype_t
enum spm_normtype_e spm_normtype_t
Norms.
spm_mtxtype_t
enum spm_mtxtype_e spm_mtxtype_t
Matrix symmetry type property.
SpmColMajor
@ SpmColMajor
Definition const.h:148
SpmDouble
@ SpmDouble
Definition const.h:64
SpmGeneral
@ SpmGeneral
Definition const.h:165
SpmCSC
@ SpmCSC
Definition const.h:73
SpmCSR
@ SpmCSR
Definition const.h:74
SpmIJV
@ SpmIJV
Definition const.h:75
SpmFrobeniusNorm
@ SpmFrobeniusNorm
Definition const.h:200
SpmInfNorm
@ SpmInfNorm
Definition const.h:201
SpmOneNorm
@ SpmOneNorm
Definition const.h:199
SpmMaxNorm
@ SpmMaxNorm
Definition const.h:202
d_spmNorm
double d_spmNorm(spm_normtype_t ntype, const spmatrix_t *spm)
Compute the norm of an spm matrix.
Definition d_spm_norm.c:1214
d_spmNormMat
double d_spmNormMat(spm_normtype_t ntype, const spmatrix_t *spm, spm_int_t n, const double *A, spm_int_t lda)
Compute the norm of a dense matrix that follows the distribution of an spm matrix.
Definition d_spm_norm.c:1278
spmatrix_s::flttype
spm_coeftype_t flttype
Definition spm.h:67
spmatrix_s::dofs
spm_int_t * dofs
Definition spm.h:85
spmatrix_s::nexp
spm_int_t nexp
Definition spm.h:78
spmatrix_s::values
void * values
Definition spm.h:92
spmatrix_s::clustnbr
int clustnbr
Definition spm.h:96
spmatrix_s::nnz
spm_int_t nnz
Definition spm.h:75
spmatrix_s::nnzexp
spm_int_t nnzexp
Definition spm.h:80
spmatrix_s::rowptr
spm_int_t * rowptr
Definition spm.h:90
spmatrix_s::fmttype
spm_fmttype_t fmttype
Definition spm.h:69
spmatrix_s::dof
spm_int_t dof
Definition spm.h:82
spmatrix_s::n
spm_int_t n
Definition spm.h:73
spmatrix_s::mtxtype
spm_mtxtype_t mtxtype
Definition spm.h:65
spmatrix_s::loc2glob
spm_int_t * loc2glob
Definition spm.h:91
spmatrix_s::layout
spm_layout_t layout
Definition spm.h:87
spmatrix_s::comm
SPM_Comm comm
Definition spm.h:97
spmatrix_s::baseval
spm_int_t baseval
Definition spm.h:71
spmatrix_s::colptr
spm_int_t * colptr
Definition spm.h:89
spmatrix_s::gNexp
spm_int_t gNexp
Definition spm.h:77
spmatrix_s::replicated
int replicated
Definition spm.h:98
spm_int_t
int spm_int_t
The main integer datatype used in spm arrays.
Definition datatypes.h:70
spmatrix_s
The sparse matrix data structure.
Definition spm.h:64