clanhe.c
Go to the documentation of this file.
00001 /* clanhe.f -- translated by f2c (version 20061008).
00002    You must link the resulting object file with libf2c:
00003         on Microsoft Windows system, link with libf2c.lib;
00004         on Linux or Unix systems, link with .../path/to/libf2c.a -lm
00005         or, if you install libf2c.a in a standard place, with -lf2c -lm
00006         -- in that order, at the end of the command line, as in
00007                 cc *.o -lf2c -lm
00008         Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
00009 
00010                 http://www.netlib.org/f2c/libf2c.zip
00011 */
00012 
00013 #include "f2c.h"
00014 #include "blaswrap.h"
00015 
00016 /* Table of constant values */
00017 
00018 static integer c__1 = 1;
00019 
00020 doublereal clanhe_(char *norm, char *uplo, integer *n, complex *a, integer *
00021         lda, real *work)
00022 {
00023     /* System generated locals */
00024     integer a_dim1, a_offset, i__1, i__2;
00025     real ret_val, r__1, r__2, r__3;
00026 
00027     /* Builtin functions */
00028     double c_abs(complex *), sqrt(doublereal);
00029 
00030     /* Local variables */
00031     integer i__, j;
00032     real sum, absa, scale;
00033     extern logical lsame_(char *, char *);
00034     real value;
00035     extern /* Subroutine */ int classq_(integer *, complex *, integer *, real 
00036             *, real *);
00037 
00038 
00039 /*  -- LAPACK auxiliary routine (version 3.2) -- */
00040 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
00041 /*     November 2006 */
00042 
00043 /*     .. Scalar Arguments .. */
00044 /*     .. */
00045 /*     .. Array Arguments .. */
00046 /*     .. */
00047 
00048 /*  Purpose */
00049 /*  ======= */
00050 
00051 /*  CLANHE  returns the value of the one norm,  or the Frobenius norm, or */
00052 /*  the  infinity norm,  or the  element of  largest absolute value  of a */
00053 /*  complex hermitian matrix A. */
00054 
00055 /*  Description */
00056 /*  =========== */
00057 
00058 /*  CLANHE returns the value */
00059 
00060 /*     CLANHE = ( max(abs(A(i,j))), NORM = 'M' or 'm' */
00061 /*              ( */
00062 /*              ( norm1(A),         NORM = '1', 'O' or 'o' */
00063 /*              ( */
00064 /*              ( normI(A),         NORM = 'I' or 'i' */
00065 /*              ( */
00066 /*              ( normF(A),         NORM = 'F', 'f', 'E' or 'e' */
00067 
00068 /*  where  norm1  denotes the  one norm of a matrix (maximum column sum), */
00069 /*  normI  denotes the  infinity norm  of a matrix  (maximum row sum) and */
00070 /*  normF  denotes the  Frobenius norm of a matrix (square root of sum of */
00071 /*  squares).  Note that  max(abs(A(i,j)))  is not a consistent matrix norm. */
00072 
00073 /*  Arguments */
00074 /*  ========= */
00075 
00076 /*  NORM    (input) CHARACTER*1 */
00077 /*          Specifies the value to be returned in CLANHE as described */
00078 /*          above. */
00079 
00080 /*  UPLO    (input) CHARACTER*1 */
00081 /*          Specifies whether the upper or lower triangular part of the */
00082 /*          hermitian matrix A is to be referenced. */
00083 /*          = 'U':  Upper triangular part of A is referenced */
00084 /*          = 'L':  Lower triangular part of A is referenced */
00085 
00086 /*  N       (input) INTEGER */
00087 /*          The order of the matrix A.  N >= 0.  When N = 0, CLANHE is */
00088 /*          set to zero. */
00089 
00090 /*  A       (input) COMPLEX array, dimension (LDA,N) */
00091 /*          The hermitian matrix A.  If UPLO = 'U', the leading n by n */
00092 /*          upper triangular part of A contains the upper triangular part */
00093 /*          of the matrix A, and the strictly lower triangular part of A */
00094 /*          is not referenced.  If UPLO = 'L', the leading n by n lower */
00095 /*          triangular part of A contains the lower triangular part of */
00096 /*          the matrix A, and the strictly upper triangular part of A is */
00097 /*          not referenced. Note that the imaginary parts of the diagonal */
00098 /*          elements need not be set and are assumed to be zero. */
00099 
00100 /*  LDA     (input) INTEGER */
00101 /*          The leading dimension of the array A.  LDA >= max(N,1). */
00102 
00103 /*  WORK    (workspace) REAL array, dimension (MAX(1,LWORK)), */
00104 /*          where LWORK >= N when NORM = 'I' or '1' or 'O'; otherwise, */
00105 /*          WORK is not referenced. */
00106 
00107 /* ===================================================================== */
00108 
00109 /*     .. Parameters .. */
00110 /*     .. */
00111 /*     .. Local Scalars .. */
00112 /*     .. */
00113 /*     .. External Functions .. */
00114 /*     .. */
00115 /*     .. External Subroutines .. */
00116 /*     .. */
00117 /*     .. Intrinsic Functions .. */
00118 /*     .. */
00119 /*     .. Executable Statements .. */
00120 
00121     /* Parameter adjustments */
00122     a_dim1 = *lda;
00123     a_offset = 1 + a_dim1;
00124     a -= a_offset;
00125     --work;
00126 
00127     /* Function Body */
00128     if (*n == 0) {
00129         value = 0.f;
00130     } else if (lsame_(norm, "M")) {
00131 
00132 /*        Find max(abs(A(i,j))). */
00133 
00134         value = 0.f;
00135         if (lsame_(uplo, "U")) {
00136             i__1 = *n;
00137             for (j = 1; j <= i__1; ++j) {
00138                 i__2 = j - 1;
00139                 for (i__ = 1; i__ <= i__2; ++i__) {
00140 /* Computing MAX */
00141                     r__1 = value, r__2 = c_abs(&a[i__ + j * a_dim1]);
00142                     value = dmax(r__1,r__2);
00143 /* L10: */
00144                 }
00145 /* Computing MAX */
00146                 i__2 = j + j * a_dim1;
00147                 r__2 = value, r__3 = (r__1 = a[i__2].r, dabs(r__1));
00148                 value = dmax(r__2,r__3);
00149 /* L20: */
00150             }
00151         } else {
00152             i__1 = *n;
00153             for (j = 1; j <= i__1; ++j) {
00154 /* Computing MAX */
00155                 i__2 = j + j * a_dim1;
00156                 r__2 = value, r__3 = (r__1 = a[i__2].r, dabs(r__1));
00157                 value = dmax(r__2,r__3);
00158                 i__2 = *n;
00159                 for (i__ = j + 1; i__ <= i__2; ++i__) {
00160 /* Computing MAX */
00161                     r__1 = value, r__2 = c_abs(&a[i__ + j * a_dim1]);
00162                     value = dmax(r__1,r__2);
00163 /* L30: */
00164                 }
00165 /* L40: */
00166             }
00167         }
00168     } else if (lsame_(norm, "I") || lsame_(norm, "O") || *(unsigned char *)norm == '1') {
00169 
00170 /*        Find normI(A) ( = norm1(A), since A is hermitian). */
00171 
00172         value = 0.f;
00173         if (lsame_(uplo, "U")) {
00174             i__1 = *n;
00175             for (j = 1; j <= i__1; ++j) {
00176                 sum = 0.f;
00177                 i__2 = j - 1;
00178                 for (i__ = 1; i__ <= i__2; ++i__) {
00179                     absa = c_abs(&a[i__ + j * a_dim1]);
00180                     sum += absa;
00181                     work[i__] += absa;
00182 /* L50: */
00183                 }
00184                 i__2 = j + j * a_dim1;
00185                 work[j] = sum + (r__1 = a[i__2].r, dabs(r__1));
00186 /* L60: */
00187             }
00188             i__1 = *n;
00189             for (i__ = 1; i__ <= i__1; ++i__) {
00190 /* Computing MAX */
00191                 r__1 = value, r__2 = work[i__];
00192                 value = dmax(r__1,r__2);
00193 /* L70: */
00194             }
00195         } else {
00196             i__1 = *n;
00197             for (i__ = 1; i__ <= i__1; ++i__) {
00198                 work[i__] = 0.f;
00199 /* L80: */
00200             }
00201             i__1 = *n;
00202             for (j = 1; j <= i__1; ++j) {
00203                 i__2 = j + j * a_dim1;
00204                 sum = work[j] + (r__1 = a[i__2].r, dabs(r__1));
00205                 i__2 = *n;
00206                 for (i__ = j + 1; i__ <= i__2; ++i__) {
00207                     absa = c_abs(&a[i__ + j * a_dim1]);
00208                     sum += absa;
00209                     work[i__] += absa;
00210 /* L90: */
00211                 }
00212                 value = dmax(value,sum);
00213 /* L100: */
00214             }
00215         }
00216     } else if (lsame_(norm, "F") || lsame_(norm, "E")) {
00217 
00218 /*        Find normF(A). */
00219 
00220         scale = 0.f;
00221         sum = 1.f;
00222         if (lsame_(uplo, "U")) {
00223             i__1 = *n;
00224             for (j = 2; j <= i__1; ++j) {
00225                 i__2 = j - 1;
00226                 classq_(&i__2, &a[j * a_dim1 + 1], &c__1, &scale, &sum);
00227 /* L110: */
00228             }
00229         } else {
00230             i__1 = *n - 1;
00231             for (j = 1; j <= i__1; ++j) {
00232                 i__2 = *n - j;
00233                 classq_(&i__2, &a[j + 1 + j * a_dim1], &c__1, &scale, &sum);
00234 /* L120: */
00235             }
00236         }
00237         sum *= 2;
00238         i__1 = *n;
00239         for (i__ = 1; i__ <= i__1; ++i__) {
00240             i__2 = i__ + i__ * a_dim1;
00241             if (a[i__2].r != 0.f) {
00242                 i__2 = i__ + i__ * a_dim1;
00243                 absa = (r__1 = a[i__2].r, dabs(r__1));
00244                 if (scale < absa) {
00245 /* Computing 2nd power */
00246                     r__1 = scale / absa;
00247                     sum = sum * (r__1 * r__1) + 1.f;
00248                     scale = absa;
00249                 } else {
00250 /* Computing 2nd power */
00251                     r__1 = absa / scale;
00252                     sum += r__1 * r__1;
00253                 }
00254             }
00255 /* L130: */
00256         }
00257         value = scale * sqrt(sum);
00258     }
00259 
00260     ret_val = value;
00261     return ret_val;
00262 
00263 /*     End of CLANHE */
00264 
00265 } /* clanhe_ */


swiftnav
Author(s):
autogenerated on Sat Jun 8 2019 18:55:30