norm.c 4.6 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230
  1. /**************************************************************************
  2. **
  3. ** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
  4. **
  5. ** Meschach Library
  6. **
  7. ** This Meschach Library is provided "as is" without any express
  8. ** or implied warranty of any kind with respect to this software.
  9. ** In particular the authors shall not be liable for any direct,
  10. ** indirect, special, incidental or consequential damages arising
  11. ** in any way from use of the software.
  12. **
  13. ** Everyone is granted permission to copy, modify and redistribute this
  14. ** Meschach Library, provided:
  15. ** 1. All copies contain this copyright notice.
  16. ** 2. All modified copies shall carry a notice stating who
  17. ** made the last modification and the date of such modification.
  18. ** 3. No charge is made for this software or works derived from it.
  19. ** This clause shall not be construed as constraining other software
  20. ** distributed on the same medium as this software, nor is a
  21. ** distribution fee considered a charge.
  22. **
  23. ***************************************************************************/
  24. /*
  25. A collection of functions for computing norms: scaled and unscaled
  26. */
  27. static char rcsid[] = "$Id: norm.c,v 1.6 1994/01/13 05:34:35 des Exp $";
  28. #include <stdio.h>
  29. #include <math.h>
  30. #include "matrix.h"
  31. /* _v_norm1 -- computes (scaled) 1-norms of vectors */
  32. #ifndef ANSI_C
  33. double _v_norm1(x,scale)
  34. VEC *x, *scale;
  35. #else
  36. double _v_norm1(const VEC *x, const VEC *scale)
  37. #endif
  38. {
  39. int i, dim;
  40. Real s, sum;
  41. if ( x == (VEC *)NULL )
  42. error(E_NULL,"_v_norm1");
  43. dim = x->dim;
  44. sum = 0.0;
  45. if ( scale == (VEC *)NULL )
  46. for ( i = 0; i < dim; i++ )
  47. sum += fabs(x->ve[i]);
  48. else if ( scale->dim < dim )
  49. error(E_SIZES,"_v_norm1");
  50. else
  51. for ( i = 0; i < dim; i++ )
  52. { s = scale->ve[i];
  53. sum += ( s== 0.0 ) ? fabs(x->ve[i]) : fabs(x->ve[i]/s);
  54. }
  55. return sum;
  56. }
  57. /* square -- returns x^2 */
  58. #ifndef ANSI_C
  59. double square(x)
  60. double x;
  61. #else
  62. double square(double x)
  63. #endif
  64. { return x*x; }
  65. /* cube -- returns x^3 */
  66. #ifndef ANSI_C
  67. double cube(x)
  68. double x;
  69. #else
  70. double cube(double x)
  71. #endif
  72. { return x*x*x; }
  73. /* _v_norm2 -- computes (scaled) 2-norm (Euclidean norm) of vectors */
  74. #ifndef ANSI_C
  75. double _v_norm2(x,scale)
  76. VEC *x, *scale;
  77. #else
  78. double _v_norm2(const VEC *x, const VEC *scale)
  79. #endif
  80. {
  81. int i, dim;
  82. Real s, sum;
  83. if ( x == (VEC *)NULL )
  84. error(E_NULL,"_v_norm2");
  85. dim = x->dim;
  86. sum = 0.0;
  87. if ( scale == (VEC *)NULL )
  88. for ( i = 0; i < dim; i++ )
  89. sum += square(x->ve[i]);
  90. else if ( scale->dim < dim )
  91. error(E_SIZES,"_v_norm2");
  92. else
  93. for ( i = 0; i < dim; i++ )
  94. { s = scale->ve[i];
  95. sum += ( s== 0.0 ) ? square(x->ve[i]) :
  96. square(x->ve[i]/s);
  97. }
  98. return sqrt(sum);
  99. }
  100. #define max(a,b) ((a) > (b) ? (a) : (b))
  101. /* _v_norm_inf -- computes (scaled) infinity-norm (supremum norm) of vectors */
  102. #ifndef ANSI_C
  103. double _v_norm_inf(x,scale)
  104. VEC *x, *scale;
  105. #else
  106. double _v_norm_inf(const VEC *x, const VEC *scale)
  107. #endif
  108. {
  109. int i, dim;
  110. Real s, maxval, tmp;
  111. if ( x == (VEC *)NULL )
  112. error(E_NULL,"_v_norm_inf");
  113. dim = x->dim;
  114. maxval = 0.0;
  115. if ( scale == (VEC *)NULL )
  116. for ( i = 0; i < dim; i++ )
  117. { tmp = fabs(x->ve[i]);
  118. maxval = max(maxval,tmp);
  119. }
  120. else if ( scale->dim < dim )
  121. error(E_SIZES,"_v_norm_inf");
  122. else
  123. for ( i = 0; i < dim; i++ )
  124. { s = scale->ve[i];
  125. tmp = ( s== 0.0 ) ? fabs(x->ve[i]) : fabs(x->ve[i]/s);
  126. maxval = max(maxval,tmp);
  127. }
  128. return maxval;
  129. }
  130. /* m_norm1 -- compute matrix 1-norm -- unscaled */
  131. #ifndef ANSI_C
  132. double m_norm1(A)
  133. MAT *A;
  134. #else
  135. double m_norm1(const MAT *A)
  136. #endif
  137. {
  138. int i, j, m, n;
  139. Real maxval, sum;
  140. if ( A == (MAT *)NULL )
  141. error(E_NULL,"m_norm1");
  142. m = A->m; n = A->n;
  143. maxval = 0.0;
  144. for ( j = 0; j < n; j++ )
  145. {
  146. sum = 0.0;
  147. for ( i = 0; i < m; i ++ )
  148. sum += fabs(A->me[i][j]);
  149. maxval = max(maxval,sum);
  150. }
  151. return maxval;
  152. }
  153. /* m_norm_inf -- compute matrix infinity-norm -- unscaled */
  154. #ifndef ANSI_C
  155. double m_norm_inf(A)
  156. MAT *A;
  157. #else
  158. double m_norm_inf(const MAT *A)
  159. #endif
  160. {
  161. int i, j, m, n;
  162. Real maxval, sum;
  163. if ( A == (MAT *)NULL )
  164. error(E_NULL,"m_norm_inf");
  165. m = A->m; n = A->n;
  166. maxval = 0.0;
  167. for ( i = 0; i < m; i++ )
  168. {
  169. sum = 0.0;
  170. for ( j = 0; j < n; j ++ )
  171. sum += fabs(A->me[i][j]);
  172. maxval = max(maxval,sum);
  173. }
  174. return maxval;
  175. }
  176. /* m_norm_frob -- compute matrix frobenius-norm -- unscaled */
  177. #ifndef ANSI_C
  178. double m_norm_frob(A)
  179. MAT *A;
  180. #else
  181. double m_norm_frob(const MAT *A)
  182. #endif
  183. {
  184. int i, j, m, n;
  185. Real sum;
  186. if ( A == (MAT *)NULL )
  187. error(E_NULL,"m_norm_frob");
  188. m = A->m; n = A->n;
  189. sum = 0.0;
  190. for ( i = 0; i < m; i++ )
  191. for ( j = 0; j < n; j ++ )
  192. sum += square(A->me[i][j]);
  193. return sqrt(sum);
  194. }