1 /* 2 Compute the Frobenius norm of a matrix. 3 */ 4 #include <stdlib.h> 5 #include <stdio.h> 6 #include <math.h> 7 #include "declarations.h" 8 Fnorm(A)9double Fnorm(A) 10 struct blockmatrix A; 11 { 12 int blk; 13 double nrm; 14 double temp; 15 16 nrm=0; 17 for (blk=1; blk<=A.nblocks; blk++) 18 { 19 switch (A.blocks[blk].blockcategory) 20 { 21 case DIAG: 22 temp=norm2(A.blocks[blk].blocksize,A.blocks[blk].data.vec+1); 23 nrm += temp*temp; 24 break; 25 case MATRIX: 26 temp=norm2(A.blocks[blk].blocksize*A.blocks[blk].blocksize, 27 A.blocks[blk].data.mat); 28 nrm += temp*temp; 29 break; 30 case PACKEDMATRIX: 31 default: 32 printf("Fnorm illegal block type \n"); 33 exit(206); 34 }; 35 }; 36 37 nrm=sqrt(nrm); 38 return(nrm); 39 } 40 41 /* 42 * The Knorm is the sum of the Fnorms of the blocks. 43 */ 44 Knorm(A)45double Knorm(A) 46 struct blockmatrix A; 47 { 48 int blk; 49 double nrm; 50 double temp; 51 52 nrm=0; 53 for (blk=1; blk<=A.nblocks; blk++) 54 { 55 switch (A.blocks[blk].blockcategory) 56 { 57 case DIAG: 58 temp=norm2(A.blocks[blk].blocksize,A.blocks[blk].data.vec+1); 59 nrm += temp; 60 break; 61 case MATRIX: 62 temp=norm2(A.blocks[blk].blocksize*A.blocks[blk].blocksize, 63 A.blocks[blk].data.mat); 64 nrm += temp; 65 break; 66 case PACKEDMATRIX: 67 default: 68 printf("Fnorm illegal block type \n"); 69 exit(206); 70 }; 71 }; 72 73 return(nrm); 74 } 75 76 mat1norm(A)77double mat1norm(A) 78 struct blockmatrix A; 79 { 80 int blk; 81 double nrm; 82 double temp; 83 84 nrm=0; 85 for (blk=1; blk<=A.nblocks; blk++) 86 { 87 switch (A.blocks[blk].blockcategory) 88 { 89 case DIAG: 90 temp=norm1(A.blocks[blk].blocksize,A.blocks[blk].data.vec+1); 91 nrm += temp; 92 break; 93 case MATRIX: 94 temp=norm1(1*A.blocks[blk].blocksize*A.blocks[blk].blocksize, 95 A.blocks[blk].data.mat); 96 nrm += temp; 97 break; 98 case PACKEDMATRIX: 99 default: 100 printf("Fnorm illegal block type \n"); 101 exit(206); 102 }; 103 }; 104 105 return(nrm); 106 } 107 108 matinfnorm(A)109double matinfnorm(A) 110 struct blockmatrix A; 111 { 112 int blk; 113 int i; 114 double nrm; 115 116 nrm=0; 117 for (blk=1; blk<=A.nblocks; blk++) 118 { 119 switch (A.blocks[blk].blockcategory) 120 { 121 case DIAG: 122 for (i=1; i<=A.blocks[blk].blocksize; i++) 123 { 124 if (fabs(A.blocks[blk].data.vec[i]) > nrm) 125 nrm=fabs(A.blocks[blk].data.vec[i]); 126 }; 127 break; 128 case MATRIX: 129 for (i=0; i<A.blocks[blk].blocksize*A.blocks[blk].blocksize; i++) 130 { 131 if (fabs(A.blocks[blk].data.vec[i]) > nrm) 132 nrm=fabs(A.blocks[blk].data.vec[i]); 133 }; 134 break; 135 case PACKEDMATRIX: 136 default: 137 printf("Fnorm illegal block type \n"); 138 exit(206); 139 }; 140 }; 141 142 return(nrm); 143 } 144 145