src/LAPACK/zlanhp.f

   1       DOUBLE PRECISION FUNCTION ZLANHP( NORM, UPLO, N, AP, WORK )
   2 *
   3 *  -- LAPACK auxiliary routine (version 3.1) --
   4 *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
   5 *     November 2006
   6 *
   7 *     .. Scalar Arguments ..
   8       CHARACTER          NORM, UPLO
   9       INTEGER            N
  10 *     ..
  11 *     .. Array Arguments ..
  12       DOUBLE PRECISION   WORK( * )
  13       COMPLEX*16         AP( * )
  14 *     ..
  15 *
  16 *  Purpose
  17 *  =======
  18 *
  19 *  ZLANHP  returns the value of the one norm,  or the Frobenius norm, or
  20 *  the  infinity norm,  or the  element of  largest absolute value  of a
  21 *  complex hermitian matrix A,  supplied in packed form.
  22 *
  23 *  Description
  24 *  ===========
  25 *
  26 *  ZLANHP returns the value
  27 *
  28 *     ZLANHP = ( max(abs(A(i,j))), NORM = 'M' or 'm'
  29 *              (
  30 *              ( norm1(A),         NORM = '1', 'O' or 'o'
  31 *              (
  32 *              ( normI(A),         NORM = 'I' or 'i'
  33 *              (
  34 *              ( normF(A),         NORM = 'F', 'f', 'E' or 'e'
  35 *
  36 *  where  norm1  denotes the  one norm of a matrix (maximum column sum),
  37 *  normI  denotes the  infinity norm  of a matrix  (maximum row sum) and
  38 *  normF  denotes the  Frobenius norm of a matrix (square root of sum of
  39 *  squares).  Note that  max(abs(A(i,j)))  is not a consistent matrix norm.
  40 *
  41 *  Arguments
  42 *  =========
  43 *
  44 *  NORM    (input) CHARACTER*1
  45 *          Specifies the value to be returned in ZLANHP as described
  46 *          above.
  47 *
  48 *  UPLO    (input) CHARACTER*1
  49 *          Specifies whether the upper or lower triangular part of the
  50 *          hermitian matrix A is supplied.
  51 *          = 'U':  Upper triangular part of A is supplied
  52 *          = 'L':  Lower triangular part of A is supplied
  53 *
  54 *  N       (input) INTEGER
  55 *          The order of the matrix A.  N >= 0.  When N = 0, ZLANHP is
  56 *          set to zero.
  57 *
  58 *  AP      (input) COMPLEX*16 array, dimension (N*(N+1)/2)
  59 *          The upper or lower triangle of the hermitian matrix A, packed
  60 *          columnwise in a linear array.  The j-th column of A is stored
  61 *          in the array AP as follows:
  62 *          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
  63 *          if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.
  64 *          Note that the  imaginary parts of the diagonal elements need
  65 *          not be set and are assumed to be zero.
  66 *
  67 *  WORK    (workspace) DOUBLE PRECISION array, dimension (MAX(1,LWORK)),
  68 *          where LWORK >= N when NORM = 'I' or '1' or 'O'; otherwise,
  69 *          WORK is not referenced.
  70 *
  71 * =====================================================================
  72 *
  73 *     .. Parameters ..
  74       DOUBLE PRECISION   ONE, ZERO
  75       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
  76 *     ..
  77 *     .. Local Scalars ..
  78       INTEGER            I, J, K
  79       DOUBLE PRECISION   ABSA, SCALE, SUM, VALUE
  80 *     ..
  81 *     .. External Functions ..
  82       LOGICAL            LSAME
  83       EXTERNAL           LSAME
  84 *     ..
  85 *     .. External Subroutines ..
  86       EXTERNAL           ZLASSQ
  87 *     ..
  88 *     .. Intrinsic Functions ..
  89       INTRINSIC          ABS, DBLE, MAX, SQRT
  90 *     ..
  91 *     .. Executable Statements ..
  92 *
  93       IF( N.EQ.0 ) THEN
  94          VALUE = ZERO
  95       ELSE IF( LSAME( NORM, 'M' ) ) THEN
  96 *
  97 *        Find max(abs(A(i,j))).
  98 *
  99          VALUE = ZERO
 100          IF( LSAME( UPLO, 'U' ) ) THEN
 101             K = 0
 102             DO 20 J = 1, N
 103                DO 10 I = K + 1, K + J - 1
 104                   VALUE = MAX( VALUE, ABS( AP( I ) ) )
 105    10          CONTINUE
 106                K = K + J
 107                VALUE = MAX( VALUE, ABS( DBLE( AP( K ) ) ) )
 108    20       CONTINUE
 109          ELSE
 110             K = 1
 111             DO 40 J = 1, N
 112                VALUE = MAX( VALUE, ABS( DBLE( AP( K ) ) ) )
 113                DO 30 I = K + 1, K + N - J
 114                   VALUE = MAX( VALUE, ABS( AP( I ) ) )
 115    30          CONTINUE
 116                K = K + N - J + 1
 117    40       CONTINUE
 118          END IF
 119       ELSE IF( ( LSAME( NORM, 'I' ) ) .OR. ( LSAME( NORM, 'O' ) ) .OR.
 120      $         ( NORM.EQ.'1' ) ) THEN
 121 *
 122 *        Find normI(A) ( = norm1(A), since A is hermitian).
 123 *
 124          VALUE = ZERO
 125          K = 1
 126          IF( LSAME( UPLO, 'U' ) ) THEN
 127             DO 60 J = 1, N
 128                SUM = ZERO
 129                DO 50 I = 1, J - 1
 130                   ABSA = ABS( AP( K ) )
 131                   SUM = SUM + ABSA
 132                   WORK( I ) = WORK( I ) + ABSA
 133                   K = K + 1
 134    50          CONTINUE
 135                WORK( J ) = SUM + ABS( DBLE( AP( K ) ) )
 136                K = K + 1
 137    60       CONTINUE
 138             DO 70 I = 1, N
 139                VALUE = MAX( VALUE, WORK( I ) )
 140    70       CONTINUE
 141          ELSE
 142             DO 80 I = 1, N
 143                WORK( I ) = ZERO
 144    80       CONTINUE
 145             DO 100 J = 1, N
 146                SUM = WORK( J ) + ABS( DBLE( AP( K ) ) )
 147                K = K + 1
 148                DO 90 I = J + 1, N
 149                   ABSA = ABS( AP( K ) )
 150                   SUM = SUM + ABSA
 151                   WORK( I ) = WORK( I ) + ABSA
 152                   K = K + 1
 153    90          CONTINUE
 154                VALUE = MAX( VALUE, SUM )
 155   100       CONTINUE
 156          END IF
 157       ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN
 158 *
 159 *        Find normF(A).
 160 *
 161          SCALE = ZERO
 162          SUM = ONE
 163          K = 2
 164          IF( LSAME( UPLO, 'U' ) ) THEN
 165             DO 110 J = 2, N
 166                CALL ZLASSQ( J-1, AP( K ), 1, SCALE, SUM )
 167                K = K + J
 168   110       CONTINUE
 169          ELSE
 170             DO 120 J = 1, N - 1
 171                CALL ZLASSQ( N-J, AP( K ), 1, SCALE, SUM )
 172                K = K + N - J + 1
 173   120       CONTINUE
 174          END IF
 175          SUM = 2*SUM
 176          K = 1
 177          DO 130 I = 1, N
 178             IF( DBLE( AP( K ) ).NE.ZERO ) THEN
 179                ABSA = ABS( DBLE( AP( K ) ) )
 180                IF( SCALE.LT.ABSA ) THEN
 181                   SUM = ONE + SUM*( SCALE / ABSA )**2
 182                   SCALE = ABSA
 183                ELSE
 184                   SUM = SUM + ( ABSA / SCALE )**2
 185                END IF
 186             END IF
 187             IF( LSAME( UPLO, 'U' ) ) THEN
 188                K = K + I + 1
 189             ELSE
 190                K = K + N - I + 1
 191             END IF
 192   130    CONTINUE
 193          VALUE = SCALE*SQRT( SUM )
 194       END IF
 195 *
 196       ZLANHP = VALUE
 197       RETURN
 198 *
 199 *     End of ZLANHP
 200 *
 201       END