Skip to content

Commit 831ac96

Browse files
Adds CRSCL
1 parent 1fafb88 commit 831ac96

File tree

3 files changed

+210
-1
lines changed

3 files changed

+210
-1
lines changed

SRC/CMakeLists.txt

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -227,7 +227,7 @@ set(CLASRC
227227
cppcon.f cppequ.f cpprfs.f cppsv.f cppsvx.f cpptrf.f cpptri.f cpptrs.f
228228
cptcon.f cpteqr.f cptrfs.f cptsv.f cptsvx.f cpttrf.f cpttrs.f cptts2.f
229229
crot.f cspcon.f cspmv.f cspr.f csprfs.f cspsv.f
230-
cspsvx.f csptrf.f csptri.f csptrs.f csrscl.f cstedc.f
230+
cspsvx.f csptrf.f csptri.f csptrs.f csrscl.f crscl.f cstedc.f
231231
cstegr.f cstein.f csteqr.f csycon.f csymv.f
232232
csyr.f csyrfs.f csysv.f csysvx.f csytf2.f csytrf.f csytri.f
233233
csytri2.f csytri2x.f csyswapr.f

SRC/cgetf2.f

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -189,6 +189,7 @@ SUBROUTINE CGETF2( M, N, A, LDA, IPIV, INFO )
189189
A( J+I, J ) = A( J+I, J ) / A( J, J )
190190
20 CONTINUE
191191
END IF
192+
! CALL CRSCL( M-J, A( J, J ), A( J+1, J ), 1 )
192193
END IF
193194
*
194195
ELSE IF( INFO.EQ.0 ) THEN

SRC/crscl.f

Lines changed: 208 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,208 @@
1+
*> \brief \b CRSCL multiplies a vector by the reciprocal of a real scalar.
2+
*
3+
* =========== DOCUMENTATION ===========
4+
*
5+
* Online html documentation available at
6+
* http://www.netlib.org/lapack/explore-html/
7+
*
8+
*> \htmlonly
9+
*> Download CRSCL + dependencies
10+
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/crscl.f">
11+
*> [TGZ]</a>
12+
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/crscl.f">
13+
*> [ZIP]</a>
14+
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/crscl.f">
15+
*> [TXT]</a>
16+
*> \endhtmlonly
17+
*
18+
* Definition:
19+
* ===========
20+
*
21+
* SUBROUTINE CRSCL( N, A, X, INCX )
22+
*
23+
* .. Scalar Arguments ..
24+
* INTEGER INCX, N
25+
* COMPLEX A
26+
* ..
27+
* .. Array Arguments ..
28+
* COMPLEX X( * )
29+
* ..
30+
*
31+
*
32+
*> \par Purpose:
33+
* =============
34+
*>
35+
*> \verbatim
36+
*>
37+
*> CRSCL multiplies an n-element complex vector x by the complex scalar
38+
*> 1/a. This is done without overflow or underflow as long as
39+
*> the final result x/a does not overflow or underflow.
40+
*> \endverbatim
41+
*
42+
* Arguments:
43+
* ==========
44+
*
45+
*> \param[in] N
46+
*> \verbatim
47+
*> N is INTEGER
48+
*> The number of components of the vector x.
49+
*> \endverbatim
50+
*>
51+
*> \param[in] A
52+
*> \verbatim
53+
*> A is COMPLEX
54+
*> The scalar a which is used to divide each component of x.
55+
*> A must not be 0, or the subroutine will divide by zero.
56+
*> \endverbatim
57+
*>
58+
*> \param[in,out] X
59+
*> \verbatim
60+
*> X is COMPLEX array, dimension
61+
*> (1+(N-1)*abs(INCX))
62+
*> The n-element vector x.
63+
*> \endverbatim
64+
*>
65+
*> \param[in] INCX
66+
*> \verbatim
67+
*> INCX is INTEGER
68+
*> The increment between successive values of the vector X.
69+
*> > 0: X(1) = X(1) and X(1+(i-1)*INCX) = x(i), 1< i<= n
70+
*> \endverbatim
71+
*
72+
* Authors:
73+
* ========
74+
*
75+
*> \author Univ. of Tennessee
76+
*> \author Univ. of California Berkeley
77+
*> \author Univ. of Colorado Denver
78+
*> \author NAG Ltd.
79+
*
80+
*> \ingroup complexOTHERauxiliary
81+
*
82+
* =====================================================================
83+
SUBROUTINE CRSCL( N, A, X, INCX )
84+
*
85+
* -- LAPACK auxiliary routine --
86+
* -- LAPACK is a software package provided by Univ. of Tennessee, --
87+
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
88+
*
89+
* .. Scalar Arguments ..
90+
INTEGER INCX, N
91+
COMPLEX A
92+
* ..
93+
* .. Array Arguments ..
94+
COMPLEX X( * )
95+
* ..
96+
*
97+
* =====================================================================
98+
*
99+
* .. Parameters ..
100+
REAL ZERO, ONE
101+
PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 )
102+
* ..
103+
* .. Local Scalars ..
104+
REAL BIGNUM, SMLNUM, HUGE, AR, AI, ABSR, ABSI, UR
105+
% , UI
106+
COMPLEX INVA
107+
* ..
108+
* .. External Functions ..
109+
REAL SLAMCH
110+
COMPLEX CLADIV
111+
EXTERNAL SLAMCH, CLADIV
112+
* ..
113+
* .. External Subroutines ..
114+
EXTERNAL CSCAL, CSSCAL
115+
* ..
116+
* .. Intrinsic Functions ..
117+
INTRINSIC ABS
118+
* ..
119+
* .. Executable Statements ..
120+
*
121+
* Quick return if possible
122+
*
123+
IF( N.LE.0 )
124+
$ RETURN
125+
*
126+
* Get machine parameters
127+
*
128+
SMLNUM = SLAMCH( 'S' )
129+
BIGNUM = ONE / SMLNUM
130+
HUGE = SLAMCH( 'O' )
131+
*
132+
* Initialize constants related to A.
133+
*
134+
AR = REAL( A )
135+
AI = AIMAG( A )
136+
ABSR = ABS( AR )
137+
ABSI = ABS( AI )
138+
*
139+
IF( ABSI.EQ.ZERO ) THEN
140+
* If alpha is real, then we can use csrscl
141+
CALL CSRSCL( N, AR, X, INCX )
142+
*
143+
ELSE IF( ABSR.EQ.ZERO ) THEN
144+
* If alpha has a zero real part, then we follow the same rules as if
145+
* alpha were real.
146+
IF( ABSI.GT.BIGNUM ) THEN
147+
INVA = CMPLX( ZERO, -BIGNUM / AI )
148+
CALL CSSCAL( N, SMLNUM, X, INCX )
149+
CALL CSCAL( N, INVA, X, INCX )
150+
ELSE IF( ABSI.LT.SMLNUM ) THEN
151+
INVA = CMPLX( ZERO, -SMLNUM / AI )
152+
CALL CSCAL( N, INVA, X, INCX )
153+
CALL CSSCAL( N, BIGNUM, X, INCX )
154+
ELSE
155+
INVA = CMPLX( ZERO, -ONE / AI )
156+
CALL CSCAL( N, INVA, X, INCX )
157+
END IF
158+
*
159+
ELSE IF( (ABSR.GE.BIGNUM).OR.(ABSI.GE.BIGNUM) ) THEN
160+
* Either real or imaginary part is too large.
161+
INVA = CLADIV( CMPLX( BIGNUM, ZERO ), A )
162+
CALL CSSCAL( N, SMLNUM, X, INCX )
163+
CALL CSCAL( N, INVA, X, INCX )
164+
*
165+
ELSE
166+
* The following numbers can be computed without NaNs and zeros.
167+
* They do not overflow simultaneously.
168+
* They are the inverse of the real and imaginary parts of 1/alpha.
169+
UR = AR + AI * ( AI / AR )
170+
UI = AI + AR * ( AR / AI )
171+
*
172+
IF( (ABS( UR ).LT.SMLNUM).OR.(ABS( UI ).LT.SMLNUM) ) THEN
173+
INVA = CMPLX( SMLNUM / UR, -SMLNUM / UI )
174+
CALL CSCAL( N, INVA, X, INCX )
175+
CALL CSSCAL( N, BIGNUM, X, INCX )
176+
ELSE IF( ABS( UR ).GT.HUGE ) THEN
177+
IF( ABSR.GE.ABSI ) THEN
178+
UR = (SMLNUM * AR) + AI * (SMLNUM * (AI / AR))
179+
ELSE
180+
UR = (SMLNUM * AR) + AI * ((SMLNUM * AI) / AR)
181+
END IF
182+
INVA = CMPLX( ONE / UR, -BIGNUM / UI )
183+
CALL CSSCAL( N, SMLNUM, X, INCX )
184+
CALL CSCAL( N, INVA, X, INCX )
185+
ELSE IF( ABS( UI ).GT.HUGE ) THEN
186+
IF( ABSI.GE.ABSR ) THEN
187+
UI = (SMLNUM * AI) + AR * (SMLNUM * (AR / AI))
188+
ELSE
189+
UI = (SMLNUM * AI) + AR * ((SMLNUM * AR) / AI)
190+
END IF
191+
INVA = CMPLX( BIGNUM / UR, -ONE / UI )
192+
CALL CSSCAL( N, SMLNUM, X, INCX )
193+
CALL CSCAL( N, INVA, X, INCX )
194+
ELSE IF( (ABS( UR ).GT.BIGNUM).OR.(ABS( UI ).GT.BIGNUM) ) THEN
195+
INVA = CMPLX( BIGNUM / UR, -BIGNUM / UI )
196+
CALL CSSCAL( N, SMLNUM, X, INCX )
197+
CALL CSCAL( N, INVA, X, INCX )
198+
ELSE
199+
INVA = CMPLX( ONE / UR, -ONE / UI )
200+
CALL CSCAL( N, INVA, X, INCX )
201+
END IF
202+
END IF
203+
*
204+
RETURN
205+
*
206+
* End of CRSCL
207+
*
208+
END

0 commit comments

Comments
 (0)