~ubuntu-branches/debian/sid/octave3.0/sid

« back to all changes in this revision

Viewing changes to libcruft/lapack/spotf2.f

  • Committer: Bazaar Package Importer
  • Author(s): Rafael Laboissiere
  • Date: 2007-12-23 16:04:15 UTC
  • Revision ID: james.westby@ubuntu.com-20071223160415-n4gk468dihy22e9v
Tags: upstream-3.0.0
ImportĀ upstreamĀ versionĀ 3.0.0

Show diffs side-by-side

added added

removed removed

Lines of Context:
 
1
      SUBROUTINE SPOTF2( UPLO, N, A, LDA, INFO )
 
2
*
 
3
*  -- LAPACK routine (version 3.1) --
 
4
*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
 
5
*     November 2006
 
6
*
 
7
*     .. Scalar Arguments ..
 
8
      CHARACTER          UPLO
 
9
      INTEGER            INFO, LDA, N
 
10
*     ..
 
11
*     .. Array Arguments ..
 
12
      REAL               A( LDA, * )
 
13
*     ..
 
14
*
 
15
*  Purpose
 
16
*  =======
 
17
*
 
18
*  SPOTF2 computes the Cholesky factorization of a real symmetric
 
19
*  positive definite matrix A.
 
20
*
 
21
*  The factorization has the form
 
22
*     A = U' * U ,  if UPLO = 'U', or
 
23
*     A = L  * L',  if UPLO = 'L',
 
24
*  where U is an upper triangular matrix and L is lower triangular.
 
25
*
 
26
*  This is the unblocked version of the algorithm, calling Level 2 BLAS.
 
27
*
 
28
*  Arguments
 
29
*  =========
 
30
*
 
31
*  UPLO    (input) CHARACTER*1
 
32
*          Specifies whether the upper or lower triangular part of the
 
33
*          symmetric matrix A is stored.
 
34
*          = 'U':  Upper triangular
 
35
*          = 'L':  Lower triangular
 
36
*
 
37
*  N       (input) INTEGER
 
38
*          The order of the matrix A.  N >= 0.
 
39
*
 
40
*  A       (input/output) REAL array, dimension (LDA,N)
 
41
*          On entry, the symmetric matrix A.  If UPLO = 'U', the leading
 
42
*          n by n upper triangular part of A contains the upper
 
43
*          triangular part of the matrix A, and the strictly lower
 
44
*          triangular part of A is not referenced.  If UPLO = 'L', the
 
45
*          leading n by n lower triangular part of A contains the lower
 
46
*          triangular part of the matrix A, and the strictly upper
 
47
*          triangular part of A is not referenced.
 
48
*
 
49
*          On exit, if INFO = 0, the factor U or L from the Cholesky
 
50
*          factorization A = U'*U  or A = L*L'.
 
51
*
 
52
*  LDA     (input) INTEGER
 
53
*          The leading dimension of the array A.  LDA >= max(1,N).
 
54
*
 
55
*  INFO    (output) INTEGER
 
56
*          = 0: successful exit
 
57
*          < 0: if INFO = -k, the k-th argument had an illegal value
 
58
*          > 0: if INFO = k, the leading minor of order k is not
 
59
*               positive definite, and the factorization could not be
 
60
*               completed.
 
61
*
 
62
*  =====================================================================
 
63
*
 
64
*     .. Parameters ..
 
65
      REAL               ONE, ZERO
 
66
      PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+0 )
 
67
*     ..
 
68
*     .. Local Scalars ..
 
69
      LOGICAL            UPPER
 
70
      INTEGER            J
 
71
      REAL               AJJ
 
72
*     ..
 
73
*     .. External Functions ..
 
74
      LOGICAL            LSAME
 
75
      REAL               SDOT
 
76
      EXTERNAL           LSAME, SDOT
 
77
*     ..
 
78
*     .. External Subroutines ..
 
79
      EXTERNAL           SGEMV, SSCAL, XERBLA
 
80
*     ..
 
81
*     .. Intrinsic Functions ..
 
82
      INTRINSIC          MAX, SQRT
 
83
*     ..
 
84
*     .. Executable Statements ..
 
85
*
 
86
*     Test the input parameters.
 
87
*
 
88
      INFO = 0
 
89
      UPPER = LSAME( UPLO, 'U' )
 
90
      IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
 
91
         INFO = -1
 
92
      ELSE IF( N.LT.0 ) THEN
 
93
         INFO = -2
 
94
      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
 
95
         INFO = -4
 
96
      END IF
 
97
      IF( INFO.NE.0 ) THEN
 
98
         CALL XERBLA( 'SPOTF2', -INFO )
 
99
         RETURN
 
100
      END IF
 
101
*
 
102
*     Quick return if possible
 
103
*
 
104
      IF( N.EQ.0 )
 
105
     $   RETURN
 
106
*
 
107
      IF( UPPER ) THEN
 
108
*
 
109
*        Compute the Cholesky factorization A = U'*U.
 
110
*
 
111
         DO 10 J = 1, N
 
112
*
 
113
*           Compute U(J,J) and test for non-positive-definiteness.
 
114
*
 
115
            AJJ = A( J, J ) - SDOT( J-1, A( 1, J ), 1, A( 1, J ), 1 )
 
116
            IF( AJJ.LE.ZERO ) THEN
 
117
               A( J, J ) = AJJ
 
118
               GO TO 30
 
119
            END IF
 
120
            AJJ = SQRT( AJJ )
 
121
            A( J, J ) = AJJ
 
122
*
 
123
*           Compute elements J+1:N of row J.
 
124
*
 
125
            IF( J.LT.N ) THEN
 
126
               CALL SGEMV( 'Transpose', J-1, N-J, -ONE, A( 1, J+1 ),
 
127
     $                     LDA, A( 1, J ), 1, ONE, A( J, J+1 ), LDA )
 
128
               CALL SSCAL( N-J, ONE / AJJ, A( J, J+1 ), LDA )
 
129
            END IF
 
130
   10    CONTINUE
 
131
      ELSE
 
132
*
 
133
*        Compute the Cholesky factorization A = L*L'.
 
134
*
 
135
         DO 20 J = 1, N
 
136
*
 
137
*           Compute L(J,J) and test for non-positive-definiteness.
 
138
*
 
139
            AJJ = A( J, J ) - SDOT( J-1, A( J, 1 ), LDA, A( J, 1 ),
 
140
     $            LDA )
 
141
            IF( AJJ.LE.ZERO ) THEN
 
142
               A( J, J ) = AJJ
 
143
               GO TO 30
 
144
            END IF
 
145
            AJJ = SQRT( AJJ )
 
146
            A( J, J ) = AJJ
 
147
*
 
148
*           Compute elements J+1:N of column J.
 
149
*
 
150
            IF( J.LT.N ) THEN
 
151
               CALL SGEMV( 'No transpose', N-J, J-1, -ONE, A( J+1, 1 ),
 
152
     $                     LDA, A( J, 1 ), LDA, ONE, A( J+1, J ), 1 )
 
153
               CALL SSCAL( N-J, ONE / AJJ, A( J+1, J ), 1 )
 
154
            END IF
 
155
   20    CONTINUE
 
156
      END IF
 
157
      GO TO 40
 
158
*
 
159
   30 CONTINUE
 
160
      INFO = J
 
161
*
 
162
   40 CONTINUE
 
163
      RETURN
 
164
*
 
165
*     End of SPOTF2
 
166
*
 
167
      END