← Back to branch summary

~ubuntu-branches/ubuntu/karmic/scilab/karmic

« back to all changes in this revision

Viewing changes to routines/lapack/dpptrf.f

  • Committer: Bazaar Package Importer
  • Author(s): Torsten Werner
  • Date: 2002-03-21 16:57:43 UTC
  • Revision ID: james.westby@ubuntu.com-20020321165743-e9mv12c1tb1plztg
Tags: upstream-2.6
Import upstream version 2.6

Show diffs side-by-side

added added

removed removed

Lines of Context:
routines/lapack/dpptrf.f
 
1
      SUBROUTINE DPPTRF( UPLO, N, AP, INFO )
 
2
*
 
3
*  -- LAPACK routine (version 2.0) --
 
4
*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
 
5
*     Courant Institute, Argonne National Lab, and Rice University
 
6
*     March 31, 1993
 
7
*
 
8
*     .. Scalar Arguments ..
 
9
      CHARACTER          UPLO
 
10
      INTEGER            INFO, N
 
11
*     ..
 
12
*     .. Array Arguments ..
 
13
      DOUBLE PRECISION   AP( * )
 
14
*     ..
 
15
*
 
16
*  Purpose
 
17
*  =======
 
18
*
 
19
*  DPPTRF computes the Cholesky factorization of a real symmetric
 
20
*  positive definite matrix A stored in packed format.
 
21
*
 
22
*  The factorization has the form
 
23
*     A = U**T * U,  if UPLO = 'U', or
 
24
*     A = L  * L**T,  if UPLO = 'L',
 
25
*  where U is an upper triangular matrix and L is lower triangular.
 
26
*
 
27
*  Arguments
 
28
*  =========
 
29
*
 
30
*  UPLO    (input) CHARACTER*1
 
31
*          = 'U':  Upper triangle of A is stored;
 
32
*          = 'L':  Lower triangle of A is stored.
 
33
*
 
34
*  N       (input) INTEGER
 
35
*          The order of the matrix A.  N >= 0.
 
36
*
 
37
*  AP      (input/output) DOUBLE PRECISION array, dimension (N*(N+1)/2)
 
38
*          On entry, the upper or lower triangle of the symmetric matrix
 
39
*          A, packed columnwise in a linear array.  The j-th column of A
 
40
*          is stored in the array AP as follows:
 
41
*          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
 
42
*          if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.
 
43
*          See below for further details.
 
44
*
 
45
*          On exit, if INFO = 0, the triangular factor U or L from the
 
46
*          Cholesky factorization A = U**T*U or A = L*L**T, in the same
 
47
*          storage format as A.
 
48
*
 
49
*  INFO    (output) INTEGER
 
50
*          = 0:  successful exit
 
51
*          < 0:  if INFO = -i, the i-th argument had an illegal value
 
52
*          > 0:  if INFO = i, the leading minor of order i is not
 
53
*                positive definite, and the factorization could not be
 
54
*                completed.
 
55
*
 
56
*  Further Details
 
57
*  ======= =======
 
58
*
 
59
*  The packed storage scheme is illustrated by the following example
 
60
*  when N = 4, UPLO = 'U':
 
61
*
 
62
*  Two-dimensional storage of the symmetric matrix A:
 
63
*
 
64
*     a11 a12 a13 a14
 
65
*         a22 a23 a24
 
66
*             a33 a34     (aij = aji)
 
67
*                 a44
 
68
*
 
69
*  Packed storage of the upper triangle of A:
 
70
*
 
71
*  AP = [ a11, a12, a22, a13, a23, a33, a14, a24, a34, a44 ]
 
72
*
 
73
*  =====================================================================
 
74
*
 
75
*     .. Parameters ..
 
76
      DOUBLE PRECISION   ONE, ZERO
 
77
      PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
 
78
*     ..
 
79
*     .. Local Scalars ..
 
80
      LOGICAL            UPPER
 
81
      INTEGER            J, JC, JJ
 
82
      DOUBLE PRECISION   AJJ
 
83
*     ..
 
84
*     .. External Functions ..
 
85
      LOGICAL            LSAME
 
86
      DOUBLE PRECISION   DDOT
 
87
      EXTERNAL           LSAME, DDOT
 
88
*     ..
 
89
*     .. External Subroutines ..
 
90
      EXTERNAL           DSCAL, DSPR, DTPSV, XERBLA
 
91
*     ..
 
92
*     .. Intrinsic Functions ..
 
93
      INTRINSIC          SQRT
 
94
*     ..
 
95
*     .. Executable Statements ..
 
96
*
 
97
*     Test the input parameters.
 
98
*
 
99
      INFO = 0
 
100
      UPPER = LSAME( UPLO, 'U' )
 
101
      IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
 
102
         INFO = -1
 
103
      ELSE IF( N.LT.0 ) THEN
 
104
         INFO = -2
 
105
      END IF
 
106
      IF( INFO.NE.0 ) THEN
 
107
         CALL XERBLA( 'DPPTRF', -INFO )
 
108
         RETURN
 
109
      END IF
 
110
*
 
111
*     Quick return if possible
 
112
*
 
113
      IF( N.EQ.0 )
 
114
     $   RETURN
 
115
*
 
116
      IF( UPPER ) THEN
 
117
*
 
118
*        Compute the Cholesky factorization A = U'*U.
 
119
*
 
120
         JJ = 0
 
121
         DO 10 J = 1, N
 
122
            JC = JJ + 1
 
123
            JJ = JJ + J
 
124
*
 
125
*           Compute elements 1:J-1 of column J.
 
126
*
 
127
            IF( J.GT.1 )
 
128
     $         CALL DTPSV( 'Upper', 'Transpose', 'Non-unit', J-1, AP,
 
129
     $                     AP( JC ), 1 )
 
130
*
 
131
*           Compute U(J,J) and test for non-positive-definiteness.
 
132
*
 
133
            AJJ = AP( JJ ) - DDOT( J-1, AP( JC ), 1, AP( JC ), 1 )
 
134
            IF( AJJ.LE.ZERO ) THEN
 
135
               AP( JJ ) = AJJ
 
136
               GO TO 30
 
137
            END IF
 
138
            AP( JJ ) = SQRT( AJJ )
 
139
   10    CONTINUE
 
140
      ELSE
 
141
*
 
142
*        Compute the Cholesky factorization A = L*L'.
 
143
*
 
144
         JJ = 1
 
145
         DO 20 J = 1, N
 
146
*
 
147
*           Compute L(J,J) and test for non-positive-definiteness.
 
148
*
 
149
            AJJ = AP( JJ )
 
150
            IF( AJJ.LE.ZERO ) THEN
 
151
               AP( JJ ) = AJJ
 
152
               GO TO 30
 
153
            END IF
 
154
            AJJ = SQRT( AJJ )
 
155
            AP( JJ ) = AJJ
 
156
*
 
157
*           Compute elements J+1:N of column J and update the trailing
 
158
*           submatrix.
 
159
*
 
160
            IF( J.LT.N ) THEN
 
161
               CALL DSCAL( N-J, ONE / AJJ, AP( JJ+1 ), 1 )
 
162
               CALL DSPR( 'Lower', N-J, -ONE, AP( JJ+1 ), 1,
 
163
     $                    AP( JJ+N-J+1 ) )
 
164
               JJ = JJ + N - J + 1
 
165
            END IF
 
166
   20    CONTINUE
 
167
      END IF
 
168
      GO TO 40
 
169
*
 
170
   30 CONTINUE
 
171
      INFO = J
 
172
*
 
173
   40 CONTINUE
 
174
      RETURN
 
175
*
 
176
*     End of DPPTRF
 
177
*
 
178
      END