← Back to branch summary

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

« back to all changes in this revision

Viewing changes to deps/openlibm/slatec/orthol.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/orthol.f
 
1
*DECK ORTHOL
 
2
      SUBROUTINE ORTHOL (A, M, N, NRDA, IFLAG, IRANK, ISCALE, DIAG,
 
3
     +   KPIVOT, SCALES, COLS, CS)
 
4
C***BEGIN PROLOGUE  ORTHOL
 
5
C***SUBSIDIARY
 
6
C***PURPOSE  Subsidiary to BVSUP
 
7
C***LIBRARY   SLATEC
 
8
C***TYPE      SINGLE PRECISION (ORTHOL-S)
 
9
C***AUTHOR  Watts, H. A., (SNLA)
 
10
C***DESCRIPTION
 
11
C
 
12
C   Reduction of the matrix A to upper triangular form by a sequence of
 
13
C   orthogonal HOUSEHOLDER transformations pre-multiplying A
 
14
C
 
15
C   Modeled after the ALGOL codes in the articles in the REFERENCES
 
16
C   section.
 
17
C
 
18
C **********************************************************************
 
19
C   INPUT
 
20
C **********************************************************************
 
21
C
 
22
C     A -- Contains the matrix to be decomposed, must be dimensioned
 
23
C           NRDA by N
 
24
C     M -- Number of rows in the matrix, M greater or equal to N
 
25
C     N -- Number of columns in the matrix, N greater or equal to 1
 
26
C     IFLAG -- Indicates the uncertainty in the matrix data
 
27
C             = 0 when the data is to be treated as exact
 
28
C             =-K when the data is assumed to be accurate to about
 
29
C                 K digits
 
30
C     ISCALE -- Scaling indicator
 
31
C               =-1 if the matrix A is to be pre-scaled by
 
32
C               columns when appropriate.
 
33
C               Otherwise no scaling will be attempted
 
34
C     NRDA -- Row dimension of A, NRDA greater or equal to M
 
35
C     DIAG,KPIVOT,COLS -- Arrays of length at least n used internally
 
36
C         ,CS,SCALES
 
37
C
 
38
C **********************************************************************
 
39
C   OUTPUT
 
40
C **********************************************************************
 
41
C
 
42
C     IFLAG - Status indicator
 
43
C            =1 for successful decomposition
 
44
C            =2 if improper input is detected
 
45
C            =3 if rank of the matrix is less than N
 
46
C     A -- Contains the reduced matrix in the strictly upper triangular
 
47
C          part and transformation information in the lower part
 
48
C     IRANK -- Contains the numerically determined matrix rank
 
49
C     DIAG -- Contains the diagonal elements of the reduced
 
50
C             triangular matrix
 
51
C     KPIVOT -- Contains the pivotal information, the column
 
52
C               interchanges performed on the original matrix are
 
53
C               recorded here.
 
54
C     SCALES -- Contains the column scaling parameters
 
55
C
 
56
C **********************************************************************
 
57
C
 
58
C***SEE ALSO  BVSUP
 
59
C***REFERENCES  G. Golub, Numerical methods for solving linear least
 
60
C                 squares problems, Numerische Mathematik 7, (1965),
 
61
C                 pp. 206-216.
 
62
C               P. Businger and G. Golub, Linear least squares
 
63
C                 solutions by Householder transformations, Numerische
 
64
C                 Mathematik  7, (1965), pp. 269-276.
 
65
C***ROUTINES CALLED  CSCALE, R1MACH, SDOT, XERMSG
 
66
C***REVISION HISTORY  (YYMMDD)
 
67
C   750601  DATE WRITTEN
 
68
C   890531  Changed all specific intrinsics to generic.  (WRB)
 
69
C   890831  Modified array declarations.  (WRB)
 
70
C   891214  Prologue converted to Version 4.0 format.  (BAB)
 
71
C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
 
72
C   900402  Added TYPE section.  (WRB)
 
73
C   910408  Updated the AUTHOR and REFERENCES sections.  (WRB)
 
74
C   920501  Reformatted the REFERENCES section.  (WRB)
 
75
C***END PROLOGUE  ORTHOL
 
76
      DIMENSION A(NRDA,*),DIAG(*),KPIVOT(*),COLS(*),CS(*),SCALES(*)
 
77
C
 
78
C **********************************************************************
 
79
C
 
80
C     MACHINE PRECISION (COMPUTER UNIT ROUNDOFF VALUE) IS DEFINED
 
81
C     BY THE FUNCTION R1MACH.
 
82
C
 
83
C***FIRST EXECUTABLE STATEMENT  ORTHOL
 
84
      URO = R1MACH(3)
 
85
C
 
86
C **********************************************************************
 
87
C
 
88
      IF (M .GE. N  .AND.  N .GE. 1  .AND.  NRDA .GE. M) GO TO 1
 
89
      IFLAG=2
 
90
      CALL XERMSG ('SLATEC', 'ORTHOL', 'INVALID INPUT PARAMETERS.', 2,
 
91
     +   1)
 
92
      RETURN
 
93
C
 
94
    1 ACC=10.*URO
 
95
      IF (IFLAG .LT. 0) ACC=MAX(ACC,10.**IFLAG)
 
96
      SRURO=SQRT(URO)
 
97
      IFLAG=1
 
98
      IRANK=N
 
99
C
 
100
C     COMPUTE NORM**2 OF JTH COLUMN AND A MATRIX NORM
 
101
C
 
102
      ANORM=0.
 
103
      DO 2 J=1,N
 
104
         KPIVOT(J)=J
 
105
         COLS(J)=SDOT(M,A(1,J),1,A(1,J),1)
 
106
         CS(J)=COLS(J)
 
107
         ANORM=ANORM+COLS(J)
 
108
    2 CONTINUE
 
109
C
 
110
C     PERFORM COLUMN SCALING ON A WHEN SPECIFIED
 
111
C
 
112
      CALL CSCALE(A,NRDA,M,N,COLS,CS,DUM,DUM,ANORM,SCALES,ISCALE,0)
 
113
C
 
114
      ANORM=SQRT(ANORM)
 
115
C
 
116
C
 
117
C     CONSTRUCTION OF UPPER TRIANGULAR MATRIX AND RECORDING OF
 
118
C     ORTHOGONAL TRANSFORMATIONS
 
119
C
 
120
C
 
121
      DO 50 K=1,N
 
122
         MK=M-K+1
 
123
         IF (K .EQ. N) GO TO 25
 
124
         KP=K+1
 
125
C
 
126
C        SEARCHING FOR PIVOTAL COLUMN
 
127
C
 
128
         DO 10 J=K,N
 
129
            IF (COLS(J) .GE. SRURO*CS(J)) GO TO 5
 
130
            COLS(J)=SDOT(MK,A(K,J),1,A(K,J),1)
 
131
            CS(J)=COLS(J)
 
132
    5       IF (J .EQ. K) GO TO 7
 
133
            IF (SIGMA .GE. 0.99*COLS(J)) GO TO 10
 
134
    7       SIGMA=COLS(J)
 
135
            JCOL=J
 
136
   10    CONTINUE
 
137
         IF (JCOL .EQ. K) GO TO 25
 
138
C
 
139
C        PERFORM COLUMN INTERCHANGE
 
140
C
 
141
         L=KPIVOT(K)
 
142
         KPIVOT(K)=KPIVOT(JCOL)
 
143
         KPIVOT(JCOL)=L
 
144
         COLS(JCOL)=COLS(K)
 
145
         COLS(K)=SIGMA
 
146
         CSS=CS(K)
 
147
         CS(K)=CS(JCOL)
 
148
         CS(JCOL)=CSS
 
149
         SC=SCALES(K)
 
150
         SCALES(K)=SCALES(JCOL)
 
151
         SCALES(JCOL)=SC
 
152
         DO 20 L=1,M
 
153
            ASAVE=A(L,K)
 
154
            A(L,K)=A(L,JCOL)
 
155
   20       A(L,JCOL)=ASAVE
 
156
C
 
157
C        CHECK RANK OF THE MATRIX
 
158
C
 
159
   25    SIG=SDOT(MK,A(K,K),1,A(K,K),1)
 
160
         DIAGK=SQRT(SIG)
 
161
         IF (DIAGK .GT. ACC*ANORM) GO TO 30
 
162
C
 
163
C        RANK DEFICIENT PROBLEM
 
164
         IFLAG=3
 
165
         IRANK=K-1
 
166
         CALL XERMSG ('SLATEC', 'ORTHOL',
 
167
     +      'RANK OF MATRIX IS LESS THAN THE NUMBER OF COLUMNS.', 1, 1)
 
168
         RETURN
 
169
C
 
170
C        CONSTRUCT AND APPLY TRANSFORMATION TO MATRIX A
 
171
C
 
172
   30    AKK=A(K,K)
 
173
         IF (AKK .GT. 0.) DIAGK=-DIAGK
 
174
         DIAG(K)=DIAGK
 
175
         A(K,K)=AKK-DIAGK
 
176
         IF (K .EQ. N) GO TO 50
 
177
         SAD=DIAGK*AKK-SIG
 
178
         DO 40 J=KP,N
 
179
            AS=SDOT(MK,A(K,K),1,A(K,J),1)/SAD
 
180
            DO 35 L=K,M
 
181
   35          A(L,J)=A(L,J)+AS*A(L,K)
 
182
   40       COLS(J)=COLS(J)-A(K,J)**2
 
183
   50 CONTINUE
 
184
C
 
185
C
 
186
      RETURN
 
187
      END