[BACK]Return to dlansf.f CVS log [TXT] [DIR] Up to [local] / rpl / lapack / lapack

Diff for /rpl/lapack/lapack/dlansf.f between versions 1.3 and 1.15

version 1.3, 2010/08/13 21:03:49 version 1.15, 2017/06/17 10:53:54
Line 1 Line 1
       DOUBLE PRECISION FUNCTION DLANSF( NORM, TRANSR, UPLO, N, A, WORK )  *> \brief \b DLANSF returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a symmetric matrix in RFP format.
   *
   *  =========== DOCUMENTATION ===========
 *  *
 *  -- LAPACK routine (version 3.2.2)                                    --  * Online html documentation available at
   *            http://www.netlib.org/lapack/explore-html/
 *  *
 *  -- Contributed by Fred Gustavson of the IBM Watson Research Center --  *> \htmlonly
 *  -- June 2010                                                       --  *> Download DLANSF + dependencies
   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlansf.f">
   *> [TGZ]</a>
   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlansf.f">
   *> [ZIP]</a>
   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlansf.f">
   *> [TXT]</a>
   *> \endhtmlonly
   *
   *  Definition:
   *  ===========
 *  *
   *       DOUBLE PRECISION FUNCTION DLANSF( NORM, TRANSR, UPLO, N, A, WORK )
   *
   *       .. Scalar Arguments ..
   *       CHARACTER          NORM, TRANSR, UPLO
   *       INTEGER            N
   *       ..
   *       .. Array Arguments ..
   *       DOUBLE PRECISION   A( 0: * ), WORK( 0: * )
   *       ..
   *
   *
   *> \par Purpose:
   *  =============
   *>
   *> \verbatim
   *>
   *> DLANSF returns the value of the one norm, or the Frobenius norm, or
   *> the infinity norm, or the element of largest absolute value of a
   *> real symmetric matrix A in RFP format.
   *> \endverbatim
   *>
   *> \return DLANSF
   *> \verbatim
   *>
   *>    DLANSF = ( max(abs(A(i,j))), NORM = 'M' or 'm'
   *>             (
   *>             ( norm1(A),         NORM = '1', 'O' or 'o'
   *>             (
   *>             ( normI(A),         NORM = 'I' or 'i'
   *>             (
   *>             ( normF(A),         NORM = 'F', 'f', 'E' or 'e'
   *>
   *> where  norm1  denotes the  one norm of a matrix (maximum column sum),
   *> normI  denotes the  infinity norm  of a matrix  (maximum row sum) and
   *> normF  denotes the  Frobenius norm of a matrix (square root of sum of
   *> squares).  Note that  max(abs(A(i,j)))  is not a  matrix norm.
   *> \endverbatim
   *
   *  Arguments:
   *  ==========
   *
   *> \param[in] NORM
   *> \verbatim
   *>          NORM is CHARACTER*1
   *>          Specifies the value to be returned in DLANSF as described
   *>          above.
   *> \endverbatim
   *>
   *> \param[in] TRANSR
   *> \verbatim
   *>          TRANSR is CHARACTER*1
   *>          Specifies whether the RFP format of A is normal or
   *>          transposed format.
   *>          = 'N':  RFP format is Normal;
   *>          = 'T':  RFP format is Transpose.
   *> \endverbatim
   *>
   *> \param[in] UPLO
   *> \verbatim
   *>          UPLO is CHARACTER*1
   *>           On entry, UPLO specifies whether the RFP matrix A came from
   *>           an upper or lower triangular matrix as follows:
   *>           = 'U': RFP A came from an upper triangular matrix;
   *>           = 'L': RFP A came from a lower triangular matrix.
   *> \endverbatim
   *>
   *> \param[in] N
   *> \verbatim
   *>          N is INTEGER
   *>          The order of the matrix A. N >= 0. When N = 0, DLANSF is
   *>          set to zero.
   *> \endverbatim
   *>
   *> \param[in] A
   *> \verbatim
   *>          A is DOUBLE PRECISION array, dimension ( N*(N+1)/2 );
   *>          On entry, the upper (if UPLO = 'U') or lower (if UPLO = 'L')
   *>          part of the symmetric matrix A stored in RFP format. See the
   *>          "Notes" below for more details.
   *>          Unchanged on exit.
   *> \endverbatim
   *>
   *> \param[out] WORK
   *> \verbatim
   *>          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK)),
   *>          where LWORK >= N when NORM = 'I' or '1' or 'O'; otherwise,
   *>          WORK is not referenced.
   *> \endverbatim
   *
   *  Authors:
   *  ========
   *
   *> \author Univ. of Tennessee
   *> \author Univ. of California Berkeley
   *> \author Univ. of Colorado Denver
   *> \author NAG Ltd.
   *
   *> \date December 2016
   *
   *> \ingroup doubleOTHERcomputational
   *
   *> \par Further Details:
   *  =====================
   *>
   *> \verbatim
   *>
   *>  We first consider Rectangular Full Packed (RFP) Format when N is
   *>  even. We give an example where N = 6.
   *>
   *>      AP is Upper             AP is Lower
   *>
   *>   00 01 02 03 04 05       00
   *>      11 12 13 14 15       10 11
   *>         22 23 24 25       20 21 22
   *>            33 34 35       30 31 32 33
   *>               44 45       40 41 42 43 44
   *>                  55       50 51 52 53 54 55
   *>
   *>
   *>  Let TRANSR = 'N'. RFP holds AP as follows:
   *>  For UPLO = 'U' the upper trapezoid A(0:5,0:2) consists of the last
   *>  three columns of AP upper. The lower triangle A(4:6,0:2) consists of
   *>  the transpose of the first three columns of AP upper.
   *>  For UPLO = 'L' the lower trapezoid A(1:6,0:2) consists of the first
   *>  three columns of AP lower. The upper triangle A(0:2,0:2) consists of
   *>  the transpose of the last three columns of AP lower.
   *>  This covers the case N even and TRANSR = 'N'.
   *>
   *>         RFP A                   RFP A
   *>
   *>        03 04 05                33 43 53
   *>        13 14 15                00 44 54
   *>        23 24 25                10 11 55
   *>        33 34 35                20 21 22
   *>        00 44 45                30 31 32
   *>        01 11 55                40 41 42
   *>        02 12 22                50 51 52
   *>
   *>  Now let TRANSR = 'T'. RFP A in both UPLO cases is just the
   *>  transpose of RFP A above. One therefore gets:
   *>
   *>
   *>           RFP A                   RFP A
   *>
   *>     03 13 23 33 00 01 02    33 00 10 20 30 40 50
   *>     04 14 24 34 44 11 12    43 44 11 21 31 41 51
   *>     05 15 25 35 45 55 22    53 54 55 22 32 42 52
   *>
   *>
   *>  We then consider Rectangular Full Packed (RFP) Format when N is
   *>  odd. We give an example where N = 5.
   *>
   *>     AP is Upper                 AP is Lower
   *>
   *>   00 01 02 03 04              00
   *>      11 12 13 14              10 11
   *>         22 23 24              20 21 22
   *>            33 34              30 31 32 33
   *>               44              40 41 42 43 44
   *>
   *>
   *>  Let TRANSR = 'N'. RFP holds AP as follows:
   *>  For UPLO = 'U' the upper trapezoid A(0:4,0:2) consists of the last
   *>  three columns of AP upper. The lower triangle A(3:4,0:1) consists of
   *>  the transpose of the first two columns of AP upper.
   *>  For UPLO = 'L' the lower trapezoid A(0:4,0:2) consists of the first
   *>  three columns of AP lower. The upper triangle A(0:1,1:2) consists of
   *>  the transpose of the last two columns of AP lower.
   *>  This covers the case N odd and TRANSR = 'N'.
   *>
   *>         RFP A                   RFP A
   *>
   *>        02 03 04                00 33 43
   *>        12 13 14                10 11 44
   *>        22 23 24                20 21 22
   *>        00 33 34                30 31 32
   *>        01 11 44                40 41 42
   *>
   *>  Now let TRANSR = 'T'. RFP A in both UPLO cases is just the
   *>  transpose of RFP A above. One therefore gets:
   *>
   *>           RFP A                   RFP A
   *>
   *>     02 12 22 00 01             00 10 20 30 40 50
   *>     03 13 23 33 11             33 11 21 31 41 51
   *>     04 14 24 34 44             43 44 22 32 42 52
   *> \endverbatim
   *
   *  =====================================================================
         DOUBLE PRECISION FUNCTION DLANSF( NORM, TRANSR, UPLO, N, A, WORK )
   *
   *  -- LAPACK computational routine (version 3.7.0) --
 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --  *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--  *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
   *     December 2016
 *  *
 *     .. Scalar Arguments ..  *     .. Scalar Arguments ..
       CHARACTER          NORM, TRANSR, UPLO        CHARACTER          NORM, TRANSR, UPLO
Line 16 Line 222
       DOUBLE PRECISION   A( 0: * ), WORK( 0: * )        DOUBLE PRECISION   A( 0: * ), WORK( 0: * )
 *     ..  *     ..
 *  *
 *  Purpose  
 *  =======  
 *  
 *  DLANSF returns the value of the one norm, or the Frobenius norm, or  
 *  the infinity norm, or the element of largest absolute value of a  
 *  real symmetric matrix A in RFP format.  
 *  
 *  Description  
 *  ===========  
 *  
 *  DLANSF returns the value  
 *  
 *     DLANSF = ( max(abs(A(i,j))), NORM = 'M' or 'm'  
 *              (  
 *              ( norm1(A),         NORM = '1', 'O' or 'o'  
 *              (  
 *              ( normI(A),         NORM = 'I' or 'i'  
 *              (  
 *              ( normF(A),         NORM = 'F', 'f', 'E' or 'e'  
 *  
 *  where  norm1  denotes the  one norm of a matrix (maximum column sum),  
 *  normI  denotes the  infinity norm  of a matrix  (maximum row sum) and  
 *  normF  denotes the  Frobenius norm of a matrix (square root of sum of  
 *  squares).  Note that  max(abs(A(i,j)))  is not a  matrix norm.  
 *  
 *  Arguments  
 *  =========  
 *  
 *  NORM    (input) CHARACTER  
 *          Specifies the value to be returned in DLANSF as described  
 *          above.  
 *  
 *  TRANSR  (input) CHARACTER  
 *          Specifies whether the RFP format of A is normal or  
 *          transposed format.  
 *          = 'N':  RFP format is Normal;  
 *          = 'T':  RFP format is Transpose.  
 *  
 *  UPLO    (input) CHARACTER  
 *           On entry, UPLO specifies whether the RFP matrix A came from  
 *           an upper or lower triangular matrix as follows:  
 *           = 'U': RFP A came from an upper triangular matrix;  
 *           = 'L': RFP A came from a lower triangular matrix.  
 *  
 *  N       (input) INTEGER  
 *          The order of the matrix A. N >= 0. When N = 0, DLANSF is  
 *          set to zero.  
 *  
 *  A       (input) DOUBLE PRECISION array, dimension ( N*(N+1)/2 );  
 *          On entry, the upper (if UPLO = 'U') or lower (if UPLO = 'L')  
 *          part of the symmetric matrix A stored in RFP format. See the  
 *          "Notes" below for more details.  
 *          Unchanged on exit.  
 *  
 *  WORK    (workspace) DOUBLE PRECISION array, dimension (MAX(1,LWORK)),  
 *          where LWORK >= N when NORM = 'I' or '1' or 'O'; otherwise,  
 *          WORK is not referenced.  
 *  
 *  Further Details  
 *  ===============  
 *  
 *  We first consider Rectangular Full Packed (RFP) Format when N is  
 *  even. We give an example where N = 6.  
 *  
 *      AP is Upper             AP is Lower  
 *  
 *   00 01 02 03 04 05       00  
 *      11 12 13 14 15       10 11  
 *         22 23 24 25       20 21 22  
 *            33 34 35       30 31 32 33  
 *               44 45       40 41 42 43 44  
 *                  55       50 51 52 53 54 55  
 *  
 *  
 *  Let TRANSR = 'N'. RFP holds AP as follows:  
 *  For UPLO = 'U' the upper trapezoid A(0:5,0:2) consists of the last  
 *  three columns of AP upper. The lower triangle A(4:6,0:2) consists of  
 *  the transpose of the first three columns of AP upper.  
 *  For UPLO = 'L' the lower trapezoid A(1:6,0:2) consists of the first  
 *  three columns of AP lower. The upper triangle A(0:2,0:2) consists of  
 *  the transpose of the last three columns of AP lower.  
 *  This covers the case N even and TRANSR = 'N'.  
 *  
 *         RFP A                   RFP A  
 *  
 *        03 04 05                33 43 53  
 *        13 14 15                00 44 54  
 *        23 24 25                10 11 55  
 *        33 34 35                20 21 22  
 *        00 44 45                30 31 32  
 *        01 11 55                40 41 42  
 *        02 12 22                50 51 52  
 *  
 *  Now let TRANSR = 'T'. RFP A in both UPLO cases is just the  
 *  transpose of RFP A above. One therefore gets:  
 *  
 *  
 *           RFP A                   RFP A  
 *  
 *     03 13 23 33 00 01 02    33 00 10 20 30 40 50  
 *     04 14 24 34 44 11 12    43 44 11 21 31 41 51  
 *     05 15 25 35 45 55 22    53 54 55 22 32 42 52  
 *  
 *  
 *  We then consider Rectangular Full Packed (RFP) Format when N is  
 *  odd. We give an example where N = 5.  
 *  
 *     AP is Upper                 AP is Lower  
 *  
 *   00 01 02 03 04              00  
 *      11 12 13 14              10 11  
 *         22 23 24              20 21 22  
 *            33 34              30 31 32 33  
 *               44              40 41 42 43 44  
 *  
 *  
 *  Let TRANSR = 'N'. RFP holds AP as follows:  
 *  For UPLO = 'U' the upper trapezoid A(0:4,0:2) consists of the last  
 *  three columns of AP upper. The lower triangle A(3:4,0:1) consists of  
 *  the transpose of the first two columns of AP upper.  
 *  For UPLO = 'L' the lower trapezoid A(0:4,0:2) consists of the first  
 *  three columns of AP lower. The upper triangle A(0:1,1:2) consists of  
 *  the transpose of the last two columns of AP lower.  
 *  This covers the case N odd and TRANSR = 'N'.  
 *  
 *         RFP A                   RFP A  
 *  
 *        02 03 04                00 33 43  
 *        12 13 14                10 11 44  
 *        22 23 24                20 21 22  
 *        00 33 34                30 31 32  
 *        01 11 44                40 41 42  
 *  
 *  Now let TRANSR = 'T'. RFP A in both UPLO cases is just the  
 *  transpose of RFP A above. One therefore gets:  
 *  
 *           RFP A                   RFP A  
 *  
 *     02 12 22 00 01             00 10 20 30 40 50  
 *     03 13 23 33 11             33 11 21 31 41 51  
 *     04 14 24 34 44             43 44 22 32 42 52  
 *  
 *  Reference  
 *  =========  
 *  
 *  =====================================================================  *  =====================================================================
 *  *
 *     .. Parameters ..  *     .. Parameters ..
Line 169 Line 230
 *     ..  *     ..
 *     .. Local Scalars ..  *     .. Local Scalars ..
       INTEGER            I, J, IFM, ILU, NOE, N1, K, L, LDA        INTEGER            I, J, IFM, ILU, NOE, N1, K, L, LDA
       DOUBLE PRECISION   SCALE, S, VALUE, AA        DOUBLE PRECISION   SCALE, S, VALUE, AA, TEMP
 *     ..  *     ..
 *     .. External Functions ..  *     .. External Functions ..
       LOGICAL            LSAME        LOGICAL            LSAME, DISNAN
       INTEGER            IDAMAX        EXTERNAL           LSAME, DISNAN
       EXTERNAL           LSAME, IDAMAX  
 *     ..  *     ..
 *     .. External Subroutines ..  *     .. External Subroutines ..
       EXTERNAL           DLASSQ        EXTERNAL           DLASSQ
Line 187 Line 247
       IF( N.EQ.0 ) THEN        IF( N.EQ.0 ) THEN
          DLANSF = ZERO           DLANSF = ZERO
          RETURN           RETURN
         ELSE IF( N.EQ.1 ) THEN
            DLANSF = ABS( A(0) )
            RETURN
       END IF        END IF
 *  *
 *     set noe = 1 if n is odd. if n is even set noe=0  *     set noe = 1 if n is odd. if n is even set noe=0
 *  *
       NOE = 1        NOE = 1
       IF( MOD( N, 2 ).EQ.0 )        IF( MOD( N, 2 ).EQ.0 )
      +   NOE = 0       $   NOE = 0
 *  *
 *     set ifm = 0 when form='T or 't' and 1 otherwise  *     set ifm = 0 when form='T or 't' and 1 otherwise
 *  *
       IFM = 1        IFM = 1
       IF( LSAME( TRANSR, 'T' ) )        IF( LSAME( TRANSR, 'T' ) )
      +   IFM = 0       $   IFM = 0
 *  *
 *     set ilu = 0 when uplo='U or 'u' and 1 otherwise  *     set ilu = 0 when uplo='U or 'u' and 1 otherwise
 *  *
       ILU = 1        ILU = 1
       IF( LSAME( UPLO, 'U' ) )        IF( LSAME( UPLO, 'U' ) )
      +   ILU = 0       $   ILU = 0
 *  *
 *     set lda = (n+1)/2 when ifm = 0  *     set lda = (n+1)/2 when ifm = 0
 *     set lda = n when ifm = 1 and noe = 1  *     set lda = n when ifm = 1 and noe = 1
Line 235 Line 298
 *           A is n by k  *           A is n by k
                DO J = 0, K - 1                 DO J = 0, K - 1
                   DO I = 0, N - 1                    DO I = 0, N - 1
                      VALUE = MAX( VALUE, ABS( A( I+J*LDA ) ) )                       TEMP = ABS( A( I+J*LDA ) )
                        IF( VALUE .LT. TEMP .OR. DISNAN( TEMP ) )
        $                    VALUE = TEMP
                   END DO                    END DO
                END DO                 END DO
             ELSE              ELSE
 *              xpose case; A is k by n  *              xpose case; A is k by n
                DO J = 0, N - 1                 DO J = 0, N - 1
                   DO I = 0, K - 1                    DO I = 0, K - 1
                      VALUE = MAX( VALUE, ABS( A( I+J*LDA ) ) )                       TEMP = ABS( A( I+J*LDA ) )
                        IF( VALUE .LT. TEMP .OR. DISNAN( TEMP ) )
        $                    VALUE = TEMP
                   END DO                    END DO
                END DO                 END DO
             END IF              END IF
Line 252 Line 319
 *              A is n+1 by k  *              A is n+1 by k
                DO J = 0, K - 1                 DO J = 0, K - 1
                   DO I = 0, N                    DO I = 0, N
                      VALUE = MAX( VALUE, ABS( A( I+J*LDA ) ) )                       TEMP = ABS( A( I+J*LDA ) )
                        IF( VALUE .LT. TEMP .OR. DISNAN( TEMP ) )
        $                    VALUE = TEMP
                   END DO                    END DO
                END DO                 END DO
             ELSE              ELSE
 *              xpose case; A is k by n+1  *              xpose case; A is k by n+1
                DO J = 0, N                 DO J = 0, N
                   DO I = 0, K - 1                    DO I = 0, K - 1
                      VALUE = MAX( VALUE, ABS( A( I+J*LDA ) ) )                       TEMP = ABS( A( I+J*LDA ) )
                        IF( VALUE .LT. TEMP .OR. DISNAN( TEMP ) )
        $                    VALUE = TEMP
                   END DO                    END DO
                END DO                 END DO
             END IF              END IF
          END IF           END IF
       ELSE IF( ( LSAME( NORM, 'I' ) ) .OR. ( LSAME( NORM, 'O' ) ) .OR.        ELSE IF( ( LSAME( NORM, 'I' ) ) .OR. ( LSAME( NORM, 'O' ) ) .OR.
      +         ( NORM.EQ.'1' ) ) THEN       $         ( NORM.EQ.'1' ) ) THEN
 *  *
 *        Find normI(A) ( = norm1(A), since A is symmetric).  *        Find normI(A) ( = norm1(A), since A is symmetric).
 *  *
Line 289 Line 360
 *                    -> A(j+k,j+k)  *                    -> A(j+k,j+k)
                      WORK( J+K ) = S + AA                       WORK( J+K ) = S + AA
                      IF( I.EQ.K+K )                       IF( I.EQ.K+K )
      +                  GO TO 10       $                  GO TO 10
                      I = I + 1                       I = I + 1
                      AA = ABS( A( I+J*LDA ) )                       AA = ABS( A( I+J*LDA ) )
 *                    -> A(j,j)  *                    -> A(j,j)
Line 305 Line 376
                      WORK( J ) = WORK( J ) + S                       WORK( J ) = WORK( J ) + S
                   END DO                    END DO
    10             CONTINUE     10             CONTINUE
                   I = IDAMAX( N, WORK, 1 )                    VALUE = WORK( 0 )
                   VALUE = WORK( I-1 )                    DO I = 1, N-1
                        TEMP = WORK( I )
                        IF( VALUE .LT. TEMP .OR. DISNAN( TEMP ) )
        $                    VALUE = TEMP
                     END DO
                ELSE                 ELSE
 *                 ilu = 1  *                 ilu = 1
                   K = K + 1                    K = K + 1
Line 343 Line 418
                      END DO                       END DO
                      WORK( J ) = WORK( J ) + S                       WORK( J ) = WORK( J ) + S
                   END DO                    END DO
                   I = IDAMAX( N, WORK, 1 )                    VALUE = WORK( 0 )
                   VALUE = WORK( I-1 )                    DO I = 1, N-1
                        TEMP = WORK( I )
                        IF( VALUE .LT. TEMP .OR. DISNAN( TEMP ) )
        $                    VALUE = TEMP
                     END DO
                END IF                 END IF
             ELSE              ELSE
 *              n is even  *              n is even
Line 377 Line 456
                      END DO                       END DO
                      WORK( J ) = WORK( J ) + S                       WORK( J ) = WORK( J ) + S
                   END DO                    END DO
                   I = IDAMAX( N, WORK, 1 )                    VALUE = WORK( 0 )
                   VALUE = WORK( I-1 )                    DO I = 1, N-1
                        TEMP = WORK( I )
                        IF( VALUE .LT. TEMP .OR. DISNAN( TEMP ) )
        $                    VALUE = TEMP
                     END DO
                ELSE                 ELSE
 *                 ilu = 1  *                 ilu = 1
                   DO I = K, N - 1                    DO I = K, N - 1
Line 411 Line 494
                      END DO                       END DO
                      WORK( J ) = WORK( J ) + S                       WORK( J ) = WORK( J ) + S
                   END DO                    END DO
                   I = IDAMAX( N, WORK, 1 )                    VALUE = WORK( 0 )
                   VALUE = WORK( I-1 )                    DO I = 1, N-1
                        TEMP = WORK( I )
                        IF( VALUE .LT. TEMP .OR. DISNAN( TEMP ) )
        $                    VALUE = TEMP
                     END DO
                END IF                 END IF
             END IF              END IF
          ELSE           ELSE
Line 473 Line 560
                      END DO                       END DO
                      WORK( J ) = WORK( J ) + S                       WORK( J ) = WORK( J ) + S
                   END DO                    END DO
                   I = IDAMAX( N, WORK, 1 )                    VALUE = WORK( 0 )
                   VALUE = WORK( I-1 )                    DO I = 1, N-1
                        TEMP = WORK( I )
                        IF( VALUE .LT. TEMP .OR. DISNAN( TEMP ) )
        $                    VALUE = TEMP
                     END DO
                ELSE                 ELSE
 *                 ilu=1  *                 ilu=1
                   K = K + 1                    K = K + 1
Line 534 Line 625
                      END DO                       END DO
                      WORK( J ) = WORK( J ) + S                       WORK( J ) = WORK( J ) + S
                   END DO                    END DO
                   I = IDAMAX( N, WORK, 1 )                    VALUE = WORK( 0 )
                   VALUE = WORK( I-1 )                    DO I = 1, N-1
                        TEMP = WORK( I )
                        IF( VALUE .LT. TEMP .OR. DISNAN( TEMP ) )
        $                    VALUE = TEMP
                     END DO
                END IF                 END IF
             ELSE              ELSE
 *              n is even  *              n is even
Line 603 Line 698
 *                 A(k-1,k-1)  *                 A(k-1,k-1)
                   S = S + AA                    S = S + AA
                   WORK( I ) = WORK( I ) + S                    WORK( I ) = WORK( I ) + S
                   I = IDAMAX( N, WORK, 1 )                    VALUE = WORK( 0 )
                   VALUE = WORK( I-1 )                    DO I = 1, N-1
                        TEMP = WORK( I )
                        IF( VALUE .LT. TEMP .OR. DISNAN( TEMP ) )
        $                    VALUE = TEMP
                     END DO
                ELSE                 ELSE
 *                 ilu=1  *                 ilu=1
                   DO I = K, N - 1                    DO I = K, N - 1
Line 672 Line 771
                      END DO                       END DO
                      WORK( J-1 ) = WORK( J-1 ) + S                       WORK( J-1 ) = WORK( J-1 ) + S
                   END DO                    END DO
                   I = IDAMAX( N, WORK, 1 )                    VALUE = WORK( 0 )
                   VALUE = WORK( I-1 )                    DO I = 1, N-1
                        TEMP = WORK( I )
                        IF( VALUE .LT. TEMP .OR. DISNAN( TEMP ) )
        $                    VALUE = TEMP
                     END DO
                END IF                 END IF
             END IF              END IF
          END IF           END IF
Line 724 Line 827
             ELSE              ELSE
 *              A is xpose  *              A is xpose
                IF( ILU.EQ.0 ) THEN                 IF( ILU.EQ.0 ) THEN
 *                 A' is upper  *                 A**T is upper
                   DO J = 1, K - 2                    DO J = 1, K - 2
                      CALL DLASSQ( J, A( 0+( K+J )*LDA ), 1, SCALE, S )                       CALL DLASSQ( J, A( 0+( K+J )*LDA ), 1, SCALE, S )
 *                    U at A(0,k)  *                    U at A(0,k)
Line 735 Line 838
                   END DO                    END DO
                   DO J = 0, K - 2                    DO J = 0, K - 2
                      CALL DLASSQ( K-J-1, A( J+1+( J+K-1 )*LDA ), 1,                       CALL DLASSQ( K-J-1, A( J+1+( J+K-1 )*LDA ), 1,
      +                            SCALE, S )       $                            SCALE, S )
 *                    L at A(0,k-1)  *                    L at A(0,k-1)
                   END DO                    END DO
                   S = S + S                    S = S + S
Line 745 Line 848
                   CALL DLASSQ( K, A( 0+( K-1 )*LDA ), LDA+1, SCALE, S )                    CALL DLASSQ( K, A( 0+( K-1 )*LDA ), LDA+1, SCALE, S )
 *                 tri L at A(0,k-1)  *                 tri L at A(0,k-1)
                ELSE                 ELSE
 *                 A' is lower  *                 A**T is lower
                   DO J = 1, K - 1                    DO J = 1, K - 1
                      CALL DLASSQ( J, A( 0+J*LDA ), 1, SCALE, S )                       CALL DLASSQ( J, A( 0+J*LDA ), 1, SCALE, S )
 *                    U at A(0,0)  *                    U at A(0,0)
Line 806 Line 909
             ELSE              ELSE
 *              A is xpose  *              A is xpose
                IF( ILU.EQ.0 ) THEN                 IF( ILU.EQ.0 ) THEN
 *                 A' is upper  *                 A**T is upper
                   DO J = 1, K - 1                    DO J = 1, K - 1
                      CALL DLASSQ( J, A( 0+( K+1+J )*LDA ), 1, SCALE, S )                       CALL DLASSQ( J, A( 0+( K+1+J )*LDA ), 1, SCALE, S )
 *                    U at A(0,k+1)  *                    U at A(0,k+1)
Line 817 Line 920
                   END DO                    END DO
                   DO J = 0, K - 2                    DO J = 0, K - 2
                      CALL DLASSQ( K-J-1, A( J+1+( J+K )*LDA ), 1, SCALE,                       CALL DLASSQ( K-J-1, A( J+1+( J+K )*LDA ), 1, SCALE,
      +                            S )       $                            S )
 *                    L at A(0,k)  *                    L at A(0,k)
                   END DO                    END DO
                   S = S + S                    S = S + S
Line 827 Line 930
                   CALL DLASSQ( K, A( 0+K*LDA ), LDA+1, SCALE, S )                    CALL DLASSQ( K, A( 0+K*LDA ), LDA+1, SCALE, S )
 *                 tri L at A(0,k)  *                 tri L at A(0,k)
                ELSE                 ELSE
 *                 A' is lower  *                 A**T is lower
                   DO J = 1, K - 1                    DO J = 1, K - 1
                      CALL DLASSQ( J, A( 0+( J+1 )*LDA ), 1, SCALE, S )                       CALL DLASSQ( J, A( 0+( J+1 )*LDA ), 1, SCALE, S )
 *                    U at A(0,1)  *                    U at A(0,1)
Removed from v.1.3  
changed lines
  Added in v.1.15

CVSweb interface <joel.bertrand@systella.fr>