← Back to branch summary

~ubuntu-branches/ubuntu/wily/julia/wily

« back to all changes in this revision

Viewing changes to deps/openlibm/slatec/du12ls.f

  • Committer: Package Import Robot
  • Author(s): Sébastien Villemot
  • Date: 2013-01-16 12:29:42 UTC
  • mfrom: (1.1.2)
  • Revision ID: package-import@ubuntu.com-20130116122942-x86e42akjq31repw
Tags: 0.0.0+20130107.gitd9656f41-1
* New upstream snashot
* No longer try to rebuild helpdb.jl.
   + debian/rules: remove helpdb.jl from build-arch rule
   + debian/control: move back python-sphinx to Build-Depends-Indep
* debian/copyright: reflect upstream changes
* Add Build-Conflicts on libatlas3-base (makes linalg tests fail)
* debian/rules: replace obsolete USE_DEBIAN makeflag by a list of
  USE_SYSTEM_* flags
* debian/rules: on non-x86 systems, use libm instead of openlibm
* dpkg-buildflags.patch: remove patch, applied upstream
* Refreshed other patches

Show diffs side-by-side

added added

removed removed

Lines of Context:
deps/openlibm/slatec/du12ls.f
 
1
*DECK DU12LS
 
2
      SUBROUTINE DU12LS (A, MDA, M, N, B, MDB, NB, MODE, KRANK, RNORM,
 
3
     +   H, W, IC, IR)
 
4
C***BEGIN PROLOGUE  DU12LS
 
5
C***SUBSIDIARY
 
6
C***PURPOSE  Subsidiary to DLLSIA
 
7
C***LIBRARY   SLATEC
 
8
C***TYPE      DOUBLE PRECISION (U12LS-S, DU12LS-D)
 
9
C***AUTHOR  (UNKNOWN)
 
10
C***DESCRIPTION
 
11
C
 
12
C        Given the Householder QR factorization of A, this
 
13
C        subroutine solves the system AX=B. If the system
 
14
C        is of reduced rank, this routine returns a solution
 
15
C        according to the selected mode.
 
16
C
 
17
C       Note - If MODE.NE.2, W is never accessed.
 
18
C
 
19
C***SEE ALSO  DLLSIA
 
20
C***ROUTINES CALLED  DAXPY, DDOT, DNRM2, DSWAP
 
21
C***REVISION HISTORY  (YYMMDD)
 
22
C   810801  DATE WRITTEN
 
23
C   890531  Changed all specific intrinsics to generic.  (WRB)
 
24
C   890831  Modified array declarations.  (WRB)
 
25
C   891214  Prologue converted to Version 4.0 format.  (BAB)
 
26
C   900328  Added TYPE section.  (WRB)
 
27
C***END PROLOGUE  DU12LS
 
28
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
 
29
      DOUBLE PRECISION DDOT,DNRM2
 
30
      DIMENSION A(MDA,*),B(MDB,*),RNORM(*),H(*),W(*)
 
31
      INTEGER IC(*),IR(*)
 
32
C***FIRST EXECUTABLE STATEMENT  DU12LS
 
33
      K=KRANK
 
34
      KP1=K+1
 
35
C
 
36
C        RANK=0
 
37
C
 
38
      IF(K.GT.0) GO TO 410
 
39
      DO 404 JB=1,NB
 
40
      RNORM(JB)=DNRM2(M,B(1,JB),1)
 
41
  404 CONTINUE
 
42
      DO 406 JB=1,NB
 
43
      DO 406 I=1,N
 
44
      B(I,JB)=0.0D0
 
45
  406 CONTINUE
 
46
      RETURN
 
47
C
 
48
C     REORDER B TO REFLECT ROW INTERCHANGES
 
49
C
 
50
  410 CONTINUE
 
51
      I=0
 
52
  412 I=I+1
 
53
      IF(I.EQ.M) GO TO 418
 
54
      J=IR(I)
 
55
      IF(J.EQ.I) GO TO 412
 
56
      IF(J.LT.0) GO TO 412
 
57
      IR(I)=-IR(I)
 
58
      DO 413 JB=1,NB
 
59
      RNORM(JB)=B(I,JB)
 
60
  413 CONTINUE
 
61
      IJ=I
 
62
  414 DO 415 JB=1,NB
 
63
      B(IJ,JB)=B(J,JB)
 
64
  415 CONTINUE
 
65
      IJ=J
 
66
      J=IR(IJ)
 
67
      IR(IJ)=-IR(IJ)
 
68
      IF(J.NE.I) GO TO 414
 
69
      DO 416 JB=1,NB
 
70
      B(IJ,JB)=RNORM(JB)
 
71
  416 CONTINUE
 
72
      GO TO 412
 
73
  418 CONTINUE
 
74
      DO 420 I=1,M
 
75
      IR(I)=ABS(IR(I))
 
76
  420 CONTINUE
 
77
C
 
78
C     APPLY HOUSEHOLDER TRANSFORMATIONS TO B
 
79
C
 
80
      DO 430 J=1,K
 
81
      TT=A(J,J)
 
82
      A(J,J)=H(J)
 
83
      DO 425 I=1,NB
 
84
      BB=-DDOT(M-J+1,A(J,J),1,B(J,I),1)/H(J)
 
85
      CALL DAXPY(M-J+1,BB,A(J,J),1,B(J,I),1)
 
86
  425 CONTINUE
 
87
      A(J,J)=TT
 
88
  430 CONTINUE
 
89
C
 
90
C        FIND NORMS OF RESIDUAL VECTOR(S)..(BEFORE OVERWRITE B)
 
91
C
 
92
      DO 440 JB=1,NB
 
93
      RNORM(JB)=DNRM2((M-K),B(KP1,JB),1)
 
94
  440 CONTINUE
 
95
C
 
96
C     BACK SOLVE UPPER TRIANGULAR R
 
97
C
 
98
      I=K
 
99
  442 DO 444 JB=1,NB
 
100
      B(I,JB)=B(I,JB)/A(I,I)
 
101
  444 CONTINUE
 
102
      IF(I.EQ.1) GO TO 450
 
103
      IM1=I-1
 
104
      DO 448 JB=1,NB
 
105
      CALL DAXPY(IM1,-B(I,JB),A(1,I),1,B(1,JB),1)
 
106
  448 CONTINUE
 
107
      I=IM1
 
108
      GO TO 442
 
109
  450 CONTINUE
 
110
C
 
111
C     RANK LT N
 
112
C
 
113
C      TRUNCATED SOLUTION
 
114
C
 
115
      IF(K.EQ.N) GO TO 480
 
116
      DO 460 JB=1,NB
 
117
      DO 460 I=KP1,N
 
118
      B(I,JB)=0.0D0
 
119
  460 CONTINUE
 
120
      IF(MODE.EQ.1) GO TO 480
 
121
C
 
122
C      MINIMAL LENGTH SOLUTION
 
123
C
 
124
      NMK=N-K
 
125
      DO 470 JB=1,NB
 
126
      DO 465 I=1,K
 
127
      TT=-DDOT(NMK,A(I,KP1),MDA,B(KP1,JB),1)/W(I)
 
128
      TT=TT-B(I,JB)
 
129
      CALL DAXPY(NMK,TT,A(I,KP1),MDA,B(KP1,JB),1)
 
130
      B(I,JB)=B(I,JB)+TT*W(I)
 
131
  465 CONTINUE
 
132
  470 CONTINUE
 
133
C
 
134
C
 
135
C     REORDER B TO REFLECT COLUMN INTERCHANGES
 
136
C
 
137
  480 CONTINUE
 
138
      I=0
 
139
  482 I=I+1
 
140
      IF(I.EQ.N) GO TO 488
 
141
      J=IC(I)
 
142
      IF(J.EQ.I) GO TO 482
 
143
      IF(J.LT.0) GO TO 482
 
144
      IC(I)=-IC(I)
 
145
  484 CALL DSWAP(NB,B(J,1),MDB,B(I,1),MDB)
 
146
      IJ=IC(J)
 
147
      IC(J)=-IC(J)
 
148
      J=IJ
 
149
      IF(J.EQ.I) GO TO 482
 
150
      GO TO 484
 
151
  488 CONTINUE
 
152
      DO 490 I=1,N
 
153
      IC(I)=ABS(IC(I))
 
154
  490 CONTINUE
 
155
C
 
156
C        SOLUTION VECTORS ARE IN FIRST N ROWS OF B(,)
 
157
C
 
158
      RETURN
 
159
      END