← Back to branch summary

~ubuntu-branches/ubuntu/maverick/openturns/maverick

« back to all changes in this revision

Viewing changes to lib/src/BlasLapack/dbdsdc.f

  • Committer: Bazaar Package Importer
  • Author(s): Fabrice Coutadeur
  • Date: 2010-05-10 17:27:55 UTC
  • mfrom: (1.1.4 upstream) (5.1.5 sid)
  • Revision ID: james.westby@ubuntu.com-20100510172755-cb5ynskknqqi5rhp
Tags: 0.13.2-2ubuntu1
* Merge with Debian testing. No changes left.
* ubuntu_fix-python-2.6.patch: fix detection of python 2.6 libs, to not use
  LOCALMODLIBS. This pulls a dependency on SSL and makes the package FTBFS.

Show diffs side-by-side

added added

removed removed

Lines of Context:
lib/src/BlasLapack/dbdsdc.f
 
1
      SUBROUTINE DBDSDC( UPLO, COMPQ, N, D, E, U, LDU, VT, LDVT, Q, IQ,
 
2
     $                   WORK, IWORK, INFO )
 
3
*
 
4
*  -- LAPACK routine (version 3.2) --
 
5
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 
6
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 
7
*     November 2006
 
8
*
 
9
*     .. Scalar Arguments ..
 
10
      CHARACTER          COMPQ, UPLO
 
11
      INTEGER            INFO, LDU, LDVT, N
 
12
*     ..
 
13
*     .. Array Arguments ..
 
14
      INTEGER            IQ( * ), IWORK( * )
 
15
      DOUBLE PRECISION   D( * ), E( * ), Q( * ), U( LDU, * ),
 
16
     $                   VT( LDVT, * ), WORK( * )
 
17
*     ..
 
18
*
 
19
*  Purpose
 
20
*  =======
 
21
*
 
22
*  DBDSDC computes the singular value decomposition (SVD) of a real
 
23
*  N-by-N (upper or lower) bidiagonal matrix B:  B = U * S * VT,
 
24
*  using a divide and conquer method, where S is a diagonal matrix
 
25
*  with non-negative diagonal elements (the singular values of B), and
 
26
*  U and VT are orthogonal matrices of left and right singular vectors,
 
27
*  respectively. DBDSDC can be used to compute all singular values,
 
28
*  and optionally, singular vectors or singular vectors in compact form.
 
29
*
 
30
*  This code makes very mild assumptions about floating point
 
31
*  arithmetic. It will work on machines with a guard digit in
 
32
*  add/subtract, or on those binary machines without guard digits
 
33
*  which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2.
 
34
*  It could conceivably fail on hexadecimal or decimal machines
 
35
*  without guard digits, but we know of none.  See DLASD3 for details.
 
36
*
 
37
*  The code currently calls DLASDQ if singular values only are desired.
 
38
*  However, it can be slightly modified to compute singular values
 
39
*  using the divide and conquer method.
 
40
*
 
41
*  Arguments
 
42
*  =========
 
43
*
 
44
*  UPLO    (input) CHARACTER*1
 
45
*          = 'U':  B is upper bidiagonal.
 
46
*          = 'L':  B is lower bidiagonal.
 
47
*
 
48
*  COMPQ   (input) CHARACTER*1
 
49
*          Specifies whether singular vectors are to be computed
 
50
*          as follows:
 
51
*          = 'N':  Compute singular values only;
 
52
*          = 'P':  Compute singular values and compute singular
 
53
*                  vectors in compact form;
 
54
*          = 'I':  Compute singular values and singular vectors.
 
55
*
 
56
*  N       (input) INTEGER
 
57
*          The order of the matrix B.  N >= 0.
 
58
*
 
59
*  D       (input/output) DOUBLE PRECISION array, dimension (N)
 
60
*          On entry, the n diagonal elements of the bidiagonal matrix B.
 
61
*          On exit, if INFO=0, the singular values of B.
 
62
*
 
63
*  E       (input/output) DOUBLE PRECISION array, dimension (N-1)
 
64
*          On entry, the elements of E contain the offdiagonal
 
65
*          elements of the bidiagonal matrix whose SVD is desired.
 
66
*          On exit, E has been destroyed.
 
67
*
 
68
*  U       (output) DOUBLE PRECISION array, dimension (LDU,N)
 
69
*          If  COMPQ = 'I', then:
 
70
*             On exit, if INFO = 0, U contains the left singular vectors
 
71
*             of the bidiagonal matrix.
 
72
*          For other values of COMPQ, U is not referenced.
 
73
*
 
74
*  LDU     (input) INTEGER
 
75
*          The leading dimension of the array U.  LDU >= 1.
 
76
*          If singular vectors are desired, then LDU >= max( 1, N ).
 
77
*
 
78
*  VT      (output) DOUBLE PRECISION array, dimension (LDVT,N)
 
79
*          If  COMPQ = 'I', then:
 
80
*             On exit, if INFO = 0, VT' contains the right singular
 
81
*             vectors of the bidiagonal matrix.
 
82
*          For other values of COMPQ, VT is not referenced.
 
83
*
 
84
*  LDVT    (input) INTEGER
 
85
*          The leading dimension of the array VT.  LDVT >= 1.
 
86
*          If singular vectors are desired, then LDVT >= max( 1, N ).
 
87
*
 
88
*  Q       (output) DOUBLE PRECISION array, dimension (LDQ)
 
89
*          If  COMPQ = 'P', then:
 
90
*             On exit, if INFO = 0, Q and IQ contain the left
 
91
*             and right singular vectors in a compact form,
 
92
*             requiring O(N log N) space instead of 2*N**2.
 
93
*             In particular, Q contains all the DOUBLE PRECISION data in
 
94
*             LDQ >= N*(11 + 2*SMLSIZ + 8*INT(LOG_2(N/(SMLSIZ+1))))
 
95
*             words of memory, where SMLSIZ is returned by ILAENV and
 
96
*             is equal to the maximum size of the subproblems at the
 
97
*             bottom of the computation tree (usually about 25).
 
98
*          For other values of COMPQ, Q is not referenced.
 
99
*
 
100
*  IQ      (output) INTEGER array, dimension (LDIQ)
 
101
*          If  COMPQ = 'P', then:
 
102
*             On exit, if INFO = 0, Q and IQ contain the left
 
103
*             and right singular vectors in a compact form,
 
104
*             requiring O(N log N) space instead of 2*N**2.
 
105
*             In particular, IQ contains all INTEGER data in
 
106
*             LDIQ >= N*(3 + 3*INT(LOG_2(N/(SMLSIZ+1))))
 
107
*             words of memory, where SMLSIZ is returned by ILAENV and
 
108
*             is equal to the maximum size of the subproblems at the
 
109
*             bottom of the computation tree (usually about 25).
 
110
*          For other values of COMPQ, IQ is not referenced.
 
111
*
 
112
*  WORK    (workspace) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
 
113
*          If COMPQ = 'N' then LWORK >= (4 * N).
 
114
*          If COMPQ = 'P' then LWORK >= (6 * N).
 
115
*          If COMPQ = 'I' then LWORK >= (3 * N**2 + 4 * N).
 
116
*
 
117
*  IWORK   (workspace) INTEGER array, dimension (8*N)
 
118
*
 
119
*  INFO    (output) INTEGER
 
120
*          = 0:  successful exit.
 
121
*          < 0:  if INFO = -i, the i-th argument had an illegal value.
 
122
*          > 0:  The algorithm failed to compute an singular value.
 
123
*                The update process of divide and conquer failed.
 
124
*
 
125
*  Further Details
 
126
*  ===============
 
127
*
 
128
*  Based on contributions by
 
129
*     Ming Gu and Huan Ren, Computer Science Division, University of
 
130
*     California at Berkeley, USA
 
131
*
 
132
*  =====================================================================
 
133
*  Changed dimension statement in comment describing E from (N) to
 
134
*  (N-1).  Sven, 17 Feb 05.
 
135
*  =====================================================================
 
136
*
 
137
*     .. Parameters ..
 
138
      DOUBLE PRECISION   ZERO, ONE, TWO
 
139
      PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0 )
 
140
*     ..
 
141
*     .. Local Scalars ..
 
142
      INTEGER            DIFL, DIFR, GIVCOL, GIVNUM, GIVPTR, I, IC,
 
143
     $                   ICOMPQ, IERR, II, IS, IU, IUPLO, IVT, J, K, KK,
 
144
     $                   MLVL, NM1, NSIZE, PERM, POLES, QSTART, SMLSIZ,
 
145
     $                   SMLSZP, SQRE, START, WSTART, Z
 
146
      DOUBLE PRECISION   CS, EPS, ORGNRM, P, R, SN
 
147
*     ..
 
148
*     .. External Functions ..
 
149
      LOGICAL            LSAME
 
150
      INTEGER            ILAENV
 
151
      DOUBLE PRECISION   DLAMCH, DLANST
 
152
      EXTERNAL           LSAME, ILAENV, DLAMCH, DLANST
 
153
*     ..
 
154
*     .. External Subroutines ..
 
155
      EXTERNAL           DCOPY, DLARTG, DLASCL, DLASD0, DLASDA, DLASDQ,
 
156
     $                   DLASET, DLASR, DSWAP, XERBLA
 
157
*     ..
 
158
*     .. Intrinsic Functions ..
 
159
      INTRINSIC          ABS, DBLE, INT, LOG, SIGN
 
160
*     ..
 
161
*     .. Executable Statements ..
 
162
*
 
163
*     Test the input parameters.
 
164
*
 
165
      INFO = 0
 
166
*
 
167
      IUPLO = 0
 
168
      IF( LSAME( UPLO, 'U' ) )
 
169
     $   IUPLO = 1
 
170
      IF( LSAME( UPLO, 'L' ) )
 
171
     $   IUPLO = 2
 
172
      IF( LSAME( COMPQ, 'N' ) ) THEN
 
173
         ICOMPQ = 0
 
174
      ELSE IF( LSAME( COMPQ, 'P' ) ) THEN
 
175
         ICOMPQ = 1
 
176
      ELSE IF( LSAME( COMPQ, 'I' ) ) THEN
 
177
         ICOMPQ = 2
 
178
      ELSE
 
179
         ICOMPQ = -1
 
180
      END IF
 
181
      IF( IUPLO.EQ.0 ) THEN
 
182
         INFO = -1
 
183
      ELSE IF( ICOMPQ.LT.0 ) THEN
 
184
         INFO = -2
 
185
      ELSE IF( N.LT.0 ) THEN
 
186
         INFO = -3
 
187
      ELSE IF( ( LDU.LT.1 ) .OR. ( ( ICOMPQ.EQ.2 ) .AND. ( LDU.LT.
 
188
     $         N ) ) ) THEN
 
189
         INFO = -7
 
190
      ELSE IF( ( LDVT.LT.1 ) .OR. ( ( ICOMPQ.EQ.2 ) .AND. ( LDVT.LT.
 
191
     $         N ) ) ) THEN
 
192
         INFO = -9
 
193
      END IF
 
194
      IF( INFO.NE.0 ) THEN
 
195
         CALL XERBLA( 'DBDSDC', -INFO )
 
196
         RETURN
 
197
      END IF
 
198
*
 
199
*     Quick return if possible
 
200
*
 
201
      IF( N.EQ.0 )
 
202
     $   RETURN
 
203
      SMLSIZ = ILAENV( 9, 'DBDSDC', ' ', 0, 0, 0, 0 )
 
204
      IF( N.EQ.1 ) THEN
 
205
         IF( ICOMPQ.EQ.1 ) THEN
 
206
            Q( 1 ) = SIGN( ONE, D( 1 ) )
 
207
            Q( 1+SMLSIZ*N ) = ONE
 
208
         ELSE IF( ICOMPQ.EQ.2 ) THEN
 
209
            U( 1, 1 ) = SIGN( ONE, D( 1 ) )
 
210
            VT( 1, 1 ) = ONE
 
211
         END IF
 
212
         D( 1 ) = ABS( D( 1 ) )
 
213
         RETURN
 
214
      END IF
 
215
      NM1 = N - 1
 
216
*
 
217
*     If matrix lower bidiagonal, rotate to be upper bidiagonal
 
218
*     by applying Givens rotations on the left
 
219
*
 
220
      WSTART = 1
 
221
      QSTART = 3
 
222
      IF( ICOMPQ.EQ.1 ) THEN
 
223
         CALL DCOPY( N, D, 1, Q( 1 ), 1 )
 
224
         CALL DCOPY( N-1, E, 1, Q( N+1 ), 1 )
 
225
      END IF
 
226
      IF( IUPLO.EQ.2 ) THEN
 
227
         QSTART = 5
 
228
         WSTART = 2*N - 1
 
229
         DO 10 I = 1, N - 1
 
230
            CALL DLARTG( D( I ), E( I ), CS, SN, R )
 
231
            D( I ) = R
 
232
            E( I ) = SN*D( I+1 )
 
233
            D( I+1 ) = CS*D( I+1 )
 
234
            IF( ICOMPQ.EQ.1 ) THEN
 
235
               Q( I+2*N ) = CS
 
236
               Q( I+3*N ) = SN
 
237
            ELSE IF( ICOMPQ.EQ.2 ) THEN
 
238
               WORK( I ) = CS
 
239
               WORK( NM1+I ) = -SN
 
240
            END IF
 
241
   10    CONTINUE
 
242
      END IF
 
243
*
 
244
*     If ICOMPQ = 0, use DLASDQ to compute the singular values.
 
245
*
 
246
      IF( ICOMPQ.EQ.0 ) THEN
 
247
         CALL DLASDQ( 'U', 0, N, 0, 0, 0, D, E, VT, LDVT, U, LDU, U,
 
248
     $                LDU, WORK( WSTART ), INFO )
 
249
         GO TO 40
 
250
      END IF
 
251
*
 
252
*     If N is smaller than the minimum divide size SMLSIZ, then solve
 
253
*     the problem with another solver.
 
254
*
 
255
      IF( N.LE.SMLSIZ ) THEN
 
256
         IF( ICOMPQ.EQ.2 ) THEN
 
257
            CALL DLASET( 'A', N, N, ZERO, ONE, U, LDU )
 
258
            CALL DLASET( 'A', N, N, ZERO, ONE, VT, LDVT )
 
259
            CALL DLASDQ( 'U', 0, N, N, N, 0, D, E, VT, LDVT, U, LDU, U,
 
260
     $                   LDU, WORK( WSTART ), INFO )
 
261
         ELSE IF( ICOMPQ.EQ.1 ) THEN
 
262
            IU = 1
 
263
            IVT = IU + N
 
264
            CALL DLASET( 'A', N, N, ZERO, ONE, Q( IU+( QSTART-1 )*N ),
 
265
     $                   N )
 
266
            CALL DLASET( 'A', N, N, ZERO, ONE, Q( IVT+( QSTART-1 )*N ),
 
267
     $                   N )
 
268
            CALL DLASDQ( 'U', 0, N, N, N, 0, D, E,
 
269
     $                   Q( IVT+( QSTART-1 )*N ), N,
 
270
     $                   Q( IU+( QSTART-1 )*N ), N,
 
271
     $                   Q( IU+( QSTART-1 )*N ), N, WORK( WSTART ),
 
272
     $                   INFO )
 
273
         END IF
 
274
         GO TO 40
 
275
      END IF
 
276
*
 
277
      IF( ICOMPQ.EQ.2 ) THEN
 
278
         CALL DLASET( 'A', N, N, ZERO, ONE, U, LDU )
 
279
         CALL DLASET( 'A', N, N, ZERO, ONE, VT, LDVT )
 
280
      END IF
 
281
*
 
282
*     Scale.
 
283
*
 
284
      ORGNRM = DLANST( 'M', N, D, E )
 
285
      IF( ORGNRM.EQ.ZERO )
 
286
     $   RETURN
 
287
      CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, N, 1, D, N, IERR )
 
288
      CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, NM1, 1, E, NM1, IERR )
 
289
*
 
290
      EPS = DLAMCH( 'Epsilon' )
 
291
*
 
292
      MLVL = INT( LOG( DBLE( N ) / DBLE( SMLSIZ+1 ) ) / LOG( TWO ) ) + 1
 
293
      SMLSZP = SMLSIZ + 1
 
294
*
 
295
      IF( ICOMPQ.EQ.1 ) THEN
 
296
         IU = 1
 
297
         IVT = 1 + SMLSIZ
 
298
         DIFL = IVT + SMLSZP
 
299
         DIFR = DIFL + MLVL
 
300
         Z = DIFR + MLVL*2
 
301
         IC = Z + MLVL
 
302
         IS = IC + 1
 
303
         POLES = IS + 1
 
304
         GIVNUM = POLES + 2*MLVL
 
305
*
 
306
         K = 1
 
307
         GIVPTR = 2
 
308
         PERM = 3
 
309
         GIVCOL = PERM + MLVL
 
310
      END IF
 
311
*
 
312
      DO 20 I = 1, N
 
313
         IF( ABS( D( I ) ).LT.EPS ) THEN
 
314
            D( I ) = SIGN( EPS, D( I ) )
 
315
         END IF
 
316
   20 CONTINUE
 
317
*
 
318
      START = 1
 
319
      SQRE = 0
 
320
*
 
321
      DO 30 I = 1, NM1
 
322
         IF( ( ABS( E( I ) ).LT.EPS ) .OR. ( I.EQ.NM1 ) ) THEN
 
323
*
 
324
*        Subproblem found. First determine its size and then
 
325
*        apply divide and conquer on it.
 
326
*
 
327
            IF( I.LT.NM1 ) THEN
 
328
*
 
329
*        A subproblem with E(I) small for I < NM1.
 
330
*
 
331
               NSIZE = I - START + 1
 
332
            ELSE IF( ABS( E( I ) ).GE.EPS ) THEN
 
333
*
 
334
*        A subproblem with E(NM1) not too small but I = NM1.
 
335
*
 
336
               NSIZE = N - START + 1
 
337
            ELSE
 
338
*
 
339
*        A subproblem with E(NM1) small. This implies an
 
340
*        1-by-1 subproblem at D(N). Solve this 1-by-1 problem
 
341
*        first.
 
342
*
 
343
               NSIZE = I - START + 1
 
344
               IF( ICOMPQ.EQ.2 ) THEN
 
345
                  U( N, N ) = SIGN( ONE, D( N ) )
 
346
                  VT( N, N ) = ONE
 
347
               ELSE IF( ICOMPQ.EQ.1 ) THEN
 
348
                  Q( N+( QSTART-1 )*N ) = SIGN( ONE, D( N ) )
 
349
                  Q( N+( SMLSIZ+QSTART-1 )*N ) = ONE
 
350
               END IF
 
351
               D( N ) = ABS( D( N ) )
 
352
            END IF
 
353
            IF( ICOMPQ.EQ.2 ) THEN
 
354
               CALL DLASD0( NSIZE, SQRE, D( START ), E( START ),
 
355
     $                      U( START, START ), LDU, VT( START, START ),
 
356
     $                      LDVT, SMLSIZ, IWORK, WORK( WSTART ), INFO )
 
357
            ELSE
 
358
               CALL DLASDA( ICOMPQ, SMLSIZ, NSIZE, SQRE, D( START ),
 
359
     $                      E( START ), Q( START+( IU+QSTART-2 )*N ), N,
 
360
     $                      Q( START+( IVT+QSTART-2 )*N ),
 
361
     $                      IQ( START+K*N ), Q( START+( DIFL+QSTART-2 )*
 
362
     $                      N ), Q( START+( DIFR+QSTART-2 )*N ),
 
363
     $                      Q( START+( Z+QSTART-2 )*N ),
 
364
     $                      Q( START+( POLES+QSTART-2 )*N ),
 
365
     $                      IQ( START+GIVPTR*N ), IQ( START+GIVCOL*N ),
 
366
     $                      N, IQ( START+PERM*N ),
 
367
     $                      Q( START+( GIVNUM+QSTART-2 )*N ),
 
368
     $                      Q( START+( IC+QSTART-2 )*N ),
 
369
     $                      Q( START+( IS+QSTART-2 )*N ),
 
370
     $                      WORK( WSTART ), IWORK, INFO )
 
371
               IF( INFO.NE.0 ) THEN
 
372
                  RETURN
 
373
               END IF
 
374
            END IF
 
375
            START = I + 1
 
376
         END IF
 
377
   30 CONTINUE
 
378
*
 
379
*     Unscale
 
380
*
 
381
      CALL DLASCL( 'G', 0, 0, ONE, ORGNRM, N, 1, D, N, IERR )
 
382
   40 CONTINUE
 
383
*
 
384
*     Use Selection Sort to minimize swaps of singular vectors
 
385
*
 
386
      DO 60 II = 2, N
 
387
         I = II - 1
 
388
         KK = I
 
389
         P = D( I )
 
390
         DO 50 J = II, N
 
391
            IF( D( J ).GT.P ) THEN
 
392
               KK = J
 
393
               P = D( J )
 
394
            END IF
 
395
   50    CONTINUE
 
396
         IF( KK.NE.I ) THEN
 
397
            D( KK ) = D( I )
 
398
            D( I ) = P
 
399
            IF( ICOMPQ.EQ.1 ) THEN
 
400
               IQ( I ) = KK
 
401
            ELSE IF( ICOMPQ.EQ.2 ) THEN
 
402
               CALL DSWAP( N, U( 1, I ), 1, U( 1, KK ), 1 )
 
403
               CALL DSWAP( N, VT( I, 1 ), LDVT, VT( KK, 1 ), LDVT )
 
404
            END IF
 
405
         ELSE IF( ICOMPQ.EQ.1 ) THEN
 
406
            IQ( I ) = I
 
407
         END IF
 
408
   60 CONTINUE
 
409
*
 
410
*     If ICOMPQ = 1, use IQ(N,1) as the indicator for UPLO
 
411
*
 
412
      IF( ICOMPQ.EQ.1 ) THEN
 
413
         IF( IUPLO.EQ.1 ) THEN
 
414
            IQ( N ) = 1
 
415
         ELSE
 
416
            IQ( N ) = 0
 
417
         END IF
 
418
      END IF
 
419
*
 
420
*     If B is lower bidiagonal, update U by those Givens rotations
 
421
*     which rotated B to be upper bidiagonal
 
422
*
 
423
      IF( ( IUPLO.EQ.2 ) .AND. ( ICOMPQ.EQ.2 ) )
 
424
     $   CALL DLASR( 'L', 'V', 'B', N, N, WORK( 1 ), WORK( N ), U, LDU )
 
425
*
 
426
      RETURN
 
427
*
 
428
*     End of DBDSDC
 
429
*
 
430
      END