zlanhp.c
Go to the documentation of this file.
00001 /* zlanhp.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 zlanhp_(char *norm, char *uplo, integer *n, doublecomplex *ap, 
00021         doublereal *work)
00022 {
00023     /* System generated locals */
00024     integer i__1, i__2;
00025     doublereal ret_val, d__1, d__2, d__3;
00026 
00027     /* Builtin functions */
00028     double z_abs(doublecomplex *), sqrt(doublereal);
00029 
00030     /* Local variables */
00031     integer i__, j, k;
00032     doublereal sum, absa, scale;
00033     extern logical lsame_(char *, char *);
00034     doublereal value;
00035     extern /* Subroutine */ int zlassq_(integer *, doublecomplex *, integer *, 
00036              doublereal *, doublereal *);
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 /*  ZLANHP  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,  supplied in packed form. */
00054 
00055 /*  Description */
00056 /*  =========== */
00057 
00058 /*  ZLANHP returns the value */
00059 
00060 /*     ZLANHP = ( 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 ZLANHP 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 supplied. */
00083 /*          = 'U':  Upper triangular part of A is supplied */
00084 /*          = 'L':  Lower triangular part of A is supplied */
00085 
00086 /*  N       (input) INTEGER */
00087 /*          The order of the matrix A.  N >= 0.  When N = 0, ZLANHP is */
00088 /*          set to zero. */
00089 
00090 /*  AP      (input) COMPLEX*16 array, dimension (N*(N+1)/2) */
00091 /*          The upper or lower triangle of the hermitian matrix A, packed */
00092 /*          columnwise in a linear array.  The j-th column of A is stored */
00093 /*          in the array AP as follows: */
00094 /*          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j; */
00095 /*          if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n. */
00096 /*          Note that the  imaginary parts of the diagonal elements need */
00097 /*          not be set and are assumed to be zero. */
00098 
00099 /*  WORK    (workspace) DOUBLE PRECISION array, dimension (MAX(1,LWORK)), */
00100 /*          where LWORK >= N when NORM = 'I' or '1' or 'O'; otherwise, */
00101 /*          WORK is not referenced. */
00102 
00103 /* ===================================================================== */
00104 
00105 /*     .. Parameters .. */
00106 /*     .. */
00107 /*     .. Local Scalars .. */
00108 /*     .. */
00109 /*     .. External Functions .. */
00110 /*     .. */
00111 /*     .. External Subroutines .. */
00112 /*     .. */
00113 /*     .. Intrinsic Functions .. */
00114 /*     .. */
00115 /*     .. Executable Statements .. */
00116 
00117     /* Parameter adjustments */
00118     --work;
00119     --ap;
00120 
00121     /* Function Body */
00122     if (*n == 0) {
00123         value = 0.;
00124     } else if (lsame_(norm, "M")) {
00125 
00126 /*        Find max(abs(A(i,j))). */
00127 
00128         value = 0.;
00129         if (lsame_(uplo, "U")) {
00130             k = 0;
00131             i__1 = *n;
00132             for (j = 1; j <= i__1; ++j) {
00133                 i__2 = k + j - 1;
00134                 for (i__ = k + 1; i__ <= i__2; ++i__) {
00135 /* Computing MAX */
00136                     d__1 = value, d__2 = z_abs(&ap[i__]);
00137                     value = max(d__1,d__2);
00138 /* L10: */
00139                 }
00140                 k += j;
00141 /* Computing MAX */
00142                 i__2 = k;
00143                 d__2 = value, d__3 = (d__1 = ap[i__2].r, abs(d__1));
00144                 value = max(d__2,d__3);
00145 /* L20: */
00146             }
00147         } else {
00148             k = 1;
00149             i__1 = *n;
00150             for (j = 1; j <= i__1; ++j) {
00151 /* Computing MAX */
00152                 i__2 = k;
00153                 d__2 = value, d__3 = (d__1 = ap[i__2].r, abs(d__1));
00154                 value = max(d__2,d__3);
00155                 i__2 = k + *n - j;
00156                 for (i__ = k + 1; i__ <= i__2; ++i__) {
00157 /* Computing MAX */
00158                     d__1 = value, d__2 = z_abs(&ap[i__]);
00159                     value = max(d__1,d__2);
00160 /* L30: */
00161                 }
00162                 k = k + *n - j + 1;
00163 /* L40: */
00164             }
00165         }
00166     } else if (lsame_(norm, "I") || lsame_(norm, "O") || *(unsigned char *)norm == '1') {
00167 
00168 /*        Find normI(A) ( = norm1(A), since A is hermitian). */
00169 
00170         value = 0.;
00171         k = 1;
00172         if (lsame_(uplo, "U")) {
00173             i__1 = *n;
00174             for (j = 1; j <= i__1; ++j) {
00175                 sum = 0.;
00176                 i__2 = j - 1;
00177                 for (i__ = 1; i__ <= i__2; ++i__) {
00178                     absa = z_abs(&ap[k]);
00179                     sum += absa;
00180                     work[i__] += absa;
00181                     ++k;
00182 /* L50: */
00183                 }
00184                 i__2 = k;
00185                 work[j] = sum + (d__1 = ap[i__2].r, abs(d__1));
00186                 ++k;
00187 /* L60: */
00188             }
00189             i__1 = *n;
00190             for (i__ = 1; i__ <= i__1; ++i__) {
00191 /* Computing MAX */
00192                 d__1 = value, d__2 = work[i__];
00193                 value = max(d__1,d__2);
00194 /* L70: */
00195             }
00196         } else {
00197             i__1 = *n;
00198             for (i__ = 1; i__ <= i__1; ++i__) {
00199                 work[i__] = 0.;
00200 /* L80: */
00201             }
00202             i__1 = *n;
00203             for (j = 1; j <= i__1; ++j) {
00204                 i__2 = k;
00205                 sum = work[j] + (d__1 = ap[i__2].r, abs(d__1));
00206                 ++k;
00207                 i__2 = *n;
00208                 for (i__ = j + 1; i__ <= i__2; ++i__) {
00209                     absa = z_abs(&ap[k]);
00210                     sum += absa;
00211                     work[i__] += absa;
00212                     ++k;
00213 /* L90: */
00214                 }
00215                 value = max(value,sum);
00216 /* L100: */
00217             }
00218         }
00219     } else if (lsame_(norm, "F") || lsame_(norm, "E")) {
00220 
00221 /*        Find normF(A). */
00222 
00223         scale = 0.;
00224         sum = 1.;
00225         k = 2;
00226         if (lsame_(uplo, "U")) {
00227             i__1 = *n;
00228             for (j = 2; j <= i__1; ++j) {
00229                 i__2 = j - 1;
00230                 zlassq_(&i__2, &ap[k], &c__1, &scale, &sum);
00231                 k += j;
00232 /* L110: */
00233             }
00234         } else {
00235             i__1 = *n - 1;
00236             for (j = 1; j <= i__1; ++j) {
00237                 i__2 = *n - j;
00238                 zlassq_(&i__2, &ap[k], &c__1, &scale, &sum);
00239                 k = k + *n - j + 1;
00240 /* L120: */
00241             }
00242         }
00243         sum *= 2;
00244         k = 1;
00245         i__1 = *n;
00246         for (i__ = 1; i__ <= i__1; ++i__) {
00247             i__2 = k;
00248             if (ap[i__2].r != 0.) {
00249                 i__2 = k;
00250                 absa = (d__1 = ap[i__2].r, abs(d__1));
00251                 if (scale < absa) {
00252 /* Computing 2nd power */
00253                     d__1 = scale / absa;
00254                     sum = sum * (d__1 * d__1) + 1.;
00255                     scale = absa;
00256                 } else {
00257 /* Computing 2nd power */
00258                     d__1 = absa / scale;
00259                     sum += d__1 * d__1;
00260                 }
00261             }
00262             if (lsame_(uplo, "U")) {
00263                 k = k + i__ + 1;
00264             } else {
00265                 k = k + *n - i__ + 1;
00266             }
00267 /* L130: */
00268         }
00269         value = scale * sqrt(sum);
00270     }
00271 
00272     ret_val = value;
00273     return ret_val;
00274 
00275 /*     End of ZLANHP */
00276 
00277 } /* zlanhp_ */


swiftnav
Author(s):
autogenerated on Sat Jun 8 2019 18:56:41