src/gmxlib/gmx_lapack/slansy.c

   1 #include <math.h>
   2
   3
   4 #include "gmx_lapack.h"
   5
   6 float
   7 F77_FUNC(slansy,SLANSY)(const char *norm, const char *uplo, int *n, float *a, int
   8         *lda, float *work)
   9 {
  10     /* System generated locals */
  11     int a_dim1, a_offset, i__1, i__2;
  12     float ret_val, d__1, d__2, d__3;
  13     int c__1 = 1;
  14
  15     /* Local variables */
  16     int i__, j;
  17     float sum, absa, scale;
  18     float value =0.0;
  19
  20     a_dim1 = *lda;
  21     a_offset = 1 + a_dim1;
  22     a -= a_offset;
  23     --work;
  24
  25     if (*n == 0) {
  26         value = 0.;
  27     } else if (*norm=='M' || *norm=='m') {
  28
  29         value = 0.;
  30         if (*uplo=='U' || *uplo=='u') {
  31             i__1 = *n;
  32             for (j = 1; j <= i__1; ++j) {
  33                 i__2 = j;
  34                 for (i__ = 1; i__ <= i__2; ++i__) {
  35                   d__2 = value;
  36                   d__3 = fabs(a[i__ + j * a_dim1]);
  37                   value = (d__2>d__3) ? d__2 : d__3;
  38                 }
  39             }
  40         } else {
  41             i__1 = *n;
  42             for (j = 1; j <= i__1; ++j) {
  43                 i__2 = *n;
  44                 for (i__ = j; i__ <= i__2; ++i__) {
  45                   d__2 = value;
  46                   d__3 = fabs(a[i__ + j * a_dim1]);
  47                     value =  (d__2>d__3) ? d__2 : d__3;
  48                 }
  49             }
  50         }
  51     } else if (*norm=='I' || *norm=='i' || *norm=='O' || *norm=='o' || *norm=='1') {
  52
  53         value = 0.;
  54         if (*uplo=='U' || *uplo=='u') {
  55             i__1 = *n;
  56             for (j = 1; j <= i__1; ++j) {
  57                 sum = 0.;
  58                 i__2 = j - 1;
  59                 for (i__ = 1; i__ <= i__2; ++i__) {
  60                     absa = fabs(a[i__ + j * a_dim1]);
  61                     sum += absa;
  62                     work[i__] += absa;
  63                 }
  64                 work[j] = sum + fabs(a[j + j * a_dim1]);
  65             }
  66             i__1 = *n;
  67             for (i__ = 1; i__ <= i__1; ++i__) {
  68                 d__1 = value, d__2 = work[i__];
  69                 value =  (d__1>d__2) ? d__1 : d__2;
  70             }
  71         } else {
  72             i__1 = *n;
  73             for (i__ = 1; i__ <= i__1; ++i__) {
  74                 work[i__] = 0.;
  75             }
  76             i__1 = *n;
  77             for (j = 1; j <= i__1; ++j) {
  78                 sum = work[j] + fabs(a[j + j * a_dim1]);
  79                 i__2 = *n;
  80                 for (i__ = j + 1; i__ <= i__2; ++i__) {
  81                     absa = fabs(a[i__ + j * a_dim1]);
  82                     sum += absa;
  83                     work[i__] += absa;
  84                 }
  85                 if(sum>value)
  86                   value = sum;
  87             }
  88         }
  89     } else if (*norm=='F' || *norm=='f' || *norm=='E' || *norm=='e') {
  90
  91         scale = 0.;
  92         sum = 1.;
  93         if (*uplo=='U' || *uplo=='u') {
  94             i__1 = *n;
  95             for (j = 2; j <= i__1; ++j) {
  96                 i__2 = j - 1;
  97                 F77_FUNC(slassq,SLASSQ)(&i__2, &a[j * a_dim1 + 1], &c__1, &scale, &sum);
  98             }
  99         } else {
 100             i__1 = *n - 1;
 101             for (j = 1; j <= i__1; ++j) {
 102                 i__2 = *n - j;
 103                 F77_FUNC(slassq,SLASSQ)(&i__2, &a[j + 1 + j * a_dim1], &c__1, &scale, &sum);
 104             }
 105         }
 106         sum *= 2;
 107         i__1 = *lda + 1;
 108         F77_FUNC(slassq,SLASSQ)(n, &a[a_offset], &i__1, &scale, &sum);
 109         value = scale * sqrt(sum);
 110     }
 111
 112     ret_val = value;
 113     return ret_val;
 114 }
 115
 116