← Back to branch summary

~ubuntu-branches/ubuntu/hoary/scilab/hoary

« back to all changes in this revision

Viewing changes to routines/lapack/zgesc2.f

  • Committer: Bazaar Package Importer
  • Author(s): Torsten Werner
  • Date: 2005-01-09 22:58:21 UTC
  • mfrom: (1.1.1 upstream)
  • Revision ID: james.westby@ubuntu.com-20050109225821-473xr8vhgugxxx5j
Tags: 3.0-12
http://bugs.debian.org/255869
changed configure.in to build scilab's own malloc.o, closes: #255869

Show diffs side-by-side

added added

removed removed

Lines of Context:
routines/lapack/zgesc2.f
 
1
      SUBROUTINE ZGESC2( N, A, LDA, RHS, IPIV, JPIV, SCALE )
 
2
*
 
3
*  -- LAPACK auxiliary routine (version 3.0) --
 
4
*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
 
5
*     Courant Institute, Argonne National Lab, and Rice University
 
6
*     June 30, 1999
 
7
*
 
8
*     .. Scalar Arguments ..
 
9
      INTEGER            LDA, N
 
10
      DOUBLE PRECISION   SCALE
 
11
*     ..
 
12
*     .. Array Arguments ..
 
13
      INTEGER            IPIV( * ), JPIV( * )
 
14
      COMPLEX*16         A( LDA, * ), RHS( * )
 
15
*     ..
 
16
*
 
17
*  Purpose
 
18
*  =======
 
19
*
 
20
*  ZGESC2 solves a system of linear equations
 
21
*
 
22
*            A * X = scale* RHS
 
23
*
 
24
*  with a general N-by-N matrix A using the LU factorization with
 
25
*  complete pivoting computed by ZGETC2.
 
26
*
 
27
*
 
28
*  Arguments
 
29
*  =========
 
30
*
 
31
*  N       (input) INTEGER
 
32
*          The number of columns of the matrix A.
 
33
*
 
34
*  A       (input) COMPLEX*16 array, dimension (LDA, N)
 
35
*          On entry, the  LU part of the factorization of the n-by-n
 
36
*          matrix A computed by ZGETC2:  A = P * L * U * Q
 
37
*
 
38
*  LDA     (input) INTEGER
 
39
*          The leading dimension of the array A.  LDA >= max(1, N).
 
40
*
 
41
*  RHS     (input/output) COMPLEX*16 array, dimension N.
 
42
*          On entry, the right hand side vector b.
 
43
*          On exit, the solution vector X.
 
44
*
 
45
*  IPIV    (iput) INTEGER array, dimension (N).
 
46
*          The pivot indices; for 1 <= i <= N, row i of the
 
47
*          matrix has been interchanged with row IPIV(i).
 
48
*
 
49
*  JPIV    (iput) INTEGER array, dimension (N).
 
50
*          The pivot indices; for 1 <= j <= N, column j of the
 
51
*          matrix has been interchanged with column JPIV(j).
 
52
*
 
53
*  SCALE    (output) DOUBLE PRECISION
 
54
*           On exit, SCALE contains the scale factor. SCALE is chosen
 
55
*           0 <= SCALE <= 1 to prevent owerflow in the solution.
 
56
*
 
57
*  Further Details
 
58
*  ===============
 
59
*
 
60
*  Based on contributions by
 
61
*     Bo Kagstrom and Peter Poromaa, Department of Computing Science,
 
62
*     Umea University, S-901 87 Umea, Sweden.
 
63
*
 
64
*  =====================================================================
 
65
*
 
66
*     .. Parameters ..
 
67
      DOUBLE PRECISION   ZERO, ONE, TWO
 
68
      PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0 )
 
69
*     ..
 
70
*     .. Local Scalars ..
 
71
      INTEGER            I, J
 
72
      DOUBLE PRECISION   BIGNUM, EPS, SMLNUM
 
73
      COMPLEX*16         TEMP
 
74
*     ..
 
75
*     .. External Subroutines ..
 
76
      EXTERNAL           ZLASWP, ZSCAL
 
77
*     ..
 
78
*     .. External Functions ..
 
79
      INTEGER            IZAMAX
 
80
      DOUBLE PRECISION   DLAMCH
 
81
      EXTERNAL           IZAMAX, DLAMCH
 
82
*     ..
 
83
*     .. Intrinsic Functions ..
 
84
      INTRINSIC          ABS, DBLE, DCMPLX
 
85
*     ..
 
86
*     .. Executable Statements ..
 
87
*
 
88
*     Set constant to control overflow
 
89
*
 
90
      EPS = DLAMCH( 'P' )
 
91
      SMLNUM = DLAMCH( 'S' ) / EPS
 
92
      BIGNUM = ONE / SMLNUM
 
93
      CALL DLABAD( SMLNUM, BIGNUM )
 
94
*
 
95
*     Apply permutations IPIV to RHS
 
96
*
 
97
      CALL ZLASWP( 1, RHS, LDA, 1, N-1, IPIV, 1 )
 
98
*
 
99
*     Solve for L part
 
100
*
 
101
      DO 20 I = 1, N - 1
 
102
         DO 10 J = I + 1, N
 
103
            RHS( J ) = RHS( J ) - A( J, I )*RHS( I )
 
104
   10    CONTINUE
 
105
   20 CONTINUE
 
106
*
 
107
*     Solve for U part
 
108
*
 
109
      SCALE = ONE
 
110
*
 
111
*     Check for scaling
 
112
*
 
113
      I = IZAMAX( N, RHS, 1 )
 
114
      IF( TWO*SMLNUM*ABS( RHS( I ) ).GT.ABS( A( N, N ) ) ) THEN
 
115
         TEMP = DCMPLX( ONE / TWO, ZERO ) / ABS( RHS( I ) )
 
116
         CALL ZSCAL( N, TEMP, RHS( 1 ), 1 )
 
117
         SCALE = SCALE*DBLE( TEMP )
 
118
      END IF
 
119
      DO 40 I = N, 1, -1
 
120
         TEMP = DCMPLX( ONE, ZERO ) / A( I, I )
 
121
         RHS( I ) = RHS( I )*TEMP
 
122
         DO 30 J = I + 1, N
 
123
            RHS( I ) = RHS( I ) - RHS( J )*( A( I, J )*TEMP )
 
124
   30    CONTINUE
 
125
   40 CONTINUE
 
126
*
 
127
*     Apply permutations JPIV to the solution (RHS)
 
128
*
 
129
      CALL ZLASWP( 1, RHS, LDA, 1, N-1, JPIV, -1 )
 
130
      RETURN
 
131
*
 
132
*     End of ZGESC2
 
133
*
 
134
      END