ScaLAPACK 2.1  2.1
ScaLAPACK: Scalable Linear Algebra PACKage
clarot.f
Go to the documentation of this file.
1  SUBROUTINE clarot( LROWS, LLEFT, LRIGHT, NL, C, S, A, LDA, XLEFT,
2  $ XRIGHT )
3 *
4 * -- LAPACK auxiliary test routine (version 3.1) --
5 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
6 * November 2006
7 *
8 * .. Scalar Arguments ..
9  LOGICAL LLEFT, LRIGHT, LROWS
10  INTEGER LDA, NL
11  COMPLEX C, S, XLEFT, XRIGHT
12 * ..
13 * .. Array Arguments ..
14  COMPLEX A( * )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * CLAROT applies a (Givens) rotation to two adjacent rows or
21 * columns, where one element of the first and/or last column/row
22 * November 2006
23 * for use on matrices stored in some format other than GE, so
24 * that elements of the matrix may be used or modified for which
25 * no array element is provided.
26 *
27 * One example is a symmetric matrix in SB format (bandwidth=4), for
28 * which UPLO='L': Two adjacent rows will have the format:
29 *
30 * row j: * * * * * . . . .
31 * row j+1: * * * * * . . . .
32 *
33 * '*' indicates elements for which storage is provided,
34 * '.' indicates elements for which no storage is provided, but
35 * are not necessarily zero; their values are determined by
36 * symmetry. ' ' indicates elements which are necessarily zero,
37 * and have no storage provided.
38 *
39 * Those columns which have two '*'s can be handled by SROT.
40 * Those columns which have no '*'s can be ignored, since as long
41 * as the Givens rotations are carefully applied to preserve
42 * symmetry, their values are determined.
43 * Those columns which have one '*' have to be handled separately,
44 * by using separate variables "p" and "q":
45 *
46 * row j: * * * * * p . . .
47 * row j+1: q * * * * * . . . .
48 *
49 * The element p would have to be set correctly, then that column
50 * is rotated, setting p to its new value. The next call to
51 * CLAROT would rotate columns j and j+1, using p, and restore
52 * symmetry. The element q would start out being zero, and be
53 * made non-zero by the rotation. Later, rotations would presumably
54 * be chosen to zero q out.
55 *
56 * Typical Calling Sequences: rotating the i-th and (i+1)-st rows.
57 * ------- ------- ---------
58 *
59 * General dense matrix:
60 *
61 * CALL CLAROT(.TRUE.,.FALSE.,.FALSE., N, C,S,
62 * A(i,1),LDA, DUMMY, DUMMY)
63 *
64 * General banded matrix in GB format:
65 *
66 * j = MAX(1, i-KL )
67 * NL = MIN( N, i+KU+1 ) + 1-j
68 * CALL CLAROT( .TRUE., i-KL.GE.1, i+KU.LT.N, NL, C,S,
69 * A(KU+i+1-j,j),LDA-1, XLEFT, XRIGHT )
70 *
71 * [ note that i+1-j is just MIN(i,KL+1) ]
72 *
73 * Symmetric banded matrix in SY format, bandwidth K,
74 * lower triangle only:
75 *
76 * j = MAX(1, i-K )
77 * NL = MIN( K+1, i ) + 1
78 * CALL CLAROT( .TRUE., i-K.GE.1, .TRUE., NL, C,S,
79 * A(i,j), LDA, XLEFT, XRIGHT )
80 *
81 * Same, but upper triangle only:
82 *
83 * NL = MIN( K+1, N-i ) + 1
84 * CALL CLAROT( .TRUE., .TRUE., i+K.LT.N, NL, C,S,
85 * A(i,i), LDA, XLEFT, XRIGHT )
86 *
87 * Symmetric banded matrix in SB format, bandwidth K,
88 * lower triangle only:
89 *
90 * [ same as for SY, except:]
91 * . . . .
92 * A(i+1-j,j), LDA-1, XLEFT, XRIGHT )
93 *
94 * [ note that i+1-j is just MIN(i,K+1) ]
95 *
96 * Same, but upper triangle only:
97 * . . .
98 * A(K+1,i), LDA-1, XLEFT, XRIGHT )
99 *
100 * Rotating columns is just the transpose of rotating rows, except
101 * for GB and SB: (rotating columns i and i+1)
102 *
103 * GB:
104 * j = MAX(1, i-KU )
105 * NL = MIN( N, i+KL+1 ) + 1-j
106 * CALL CLAROT( .TRUE., i-KU.GE.1, i+KL.LT.N, NL, C,S,
107 * A(KU+j+1-i,i),LDA-1, XTOP, XBOTTM )
108 *
109 * [note that KU+j+1-i is just MAX(1,KU+2-i)]
110 *
111 * SB: (upper triangle)
112 *
113 * . . . . . .
114 * A(K+j+1-i,i),LDA-1, XTOP, XBOTTM )
115 *
116 * SB: (lower triangle)
117 *
118 * . . . . . .
119 * A(1,i),LDA-1, XTOP, XBOTTM )
120 *
121 * Arguments
122 * =========
123 *
124 * LROWS - LOGICAL
125 * If .TRUE., then CLAROT will rotate two rows. If .FALSE.,
126 * then it will rotate two columns.
127 * Not modified.
128 *
129 * LLEFT - LOGICAL
130 * If .TRUE., then XLEFT will be used instead of the
131 * corresponding element of A for the first element in the
132 * second row (if LROWS=.FALSE.) or column (if LROWS=.TRUE.)
133 * If .FALSE., then the corresponding element of A will be
134 * used.
135 * Not modified.
136 *
137 * LRIGHT - LOGICAL
138 * If .TRUE., then XRIGHT will be used instead of the
139 * corresponding element of A for the last element in the
140 * first row (if LROWS=.FALSE.) or column (if LROWS=.TRUE.) If
141 * .FALSE., then the corresponding element of A will be used.
142 * Not modified.
143 *
144 * NL - INTEGER
145 * The length of the rows (if LROWS=.TRUE.) or columns (if
146 * LROWS=.FALSE.) to be rotated. If XLEFT and/or XRIGHT are
147 * used, the columns/rows they are in should be included in
148 * NL, e.g., if LLEFT = LRIGHT = .TRUE., then NL must be at
149 * least 2. The number of rows/columns to be rotated
150 * exclusive of those involving XLEFT and/or XRIGHT may
151 * not be negative, i.e., NL minus how many of LLEFT and
152 * LRIGHT are .TRUE. must be at least zero; if not, XERBLA
153 * will be called.
154 * Not modified.
155 *
156 * C, S - COMPLEX
157 * Specify the Givens rotation to be applied. If LROWS is
158 * true, then the matrix ( c s )
159 * ( _ _ )
160 * (-s c ) is applied from the left;
161 * if false, then the transpose (not conjugated) thereof is
162 * applied from the right. Note that in contrast to the
163 * output of CROTG or to most versions of CROT, both C and S
164 * are complex. For a Givens rotation, |C|**2 + |S|**2 should
165 * be 1, but this is not checked.
166 * Not modified.
167 *
168 * A - COMPLEX array.
169 * The array containing the rows/columns to be rotated. The
170 * first element of A should be the upper left element to
171 * be rotated.
172 * Read and modified.
173 *
174 * LDA - INTEGER
175 * The "effective" leading dimension of A. If A contains
176 * a matrix stored in GE, HE, or SY format, then this is just
177 * the leading dimension of A as dimensioned in the calling
178 * routine. If A contains a matrix stored in band (GB, HB, or
179 * SB) format, then this should be *one less* than the leading
180 * dimension used in the calling routine. Thus, if A were
181 * dimensioned A(LDA,*) in CLAROT, then A(1,j) would be the
182 * j-th element in the first of the two rows to be rotated,
183 * and A(2,j) would be the j-th in the second, regardless of
184 * how the array may be stored in the calling routine. [A
185 * cannot, however, actually be dimensioned thus, since for
186 * band format, the row number may exceed LDA, which is not
187 * legal FORTRAN.]
188 * If LROWS=.TRUE., then LDA must be at least 1, otherwise
189 * it must be at least NL minus the number of .TRUE. values
190 * in XLEFT and XRIGHT.
191 * Not modified.
192 *
193 * XLEFT - COMPLEX
194 * If LLEFT is .TRUE., then XLEFT will be used and modified
195 * instead of A(2,1) (if LROWS=.TRUE.) or A(1,2)
196 * (if LROWS=.FALSE.).
197 * Read and modified.
198 *
199 * XRIGHT - COMPLEX
200 * If LRIGHT is .TRUE., then XRIGHT will be used and modified
201 * instead of A(1,NL) (if LROWS=.TRUE.) or A(NL,1)
202 * (if LROWS=.FALSE.).
203 * Read and modified.
204 *
205 * =====================================================================
206 *
207 * .. Local Scalars ..
208  INTEGER IINC, INEXT, IX, IY, IYT, J, NT
209  COMPLEX TEMPX
210 * ..
211 * .. Local Arrays ..
212  COMPLEX XT( 2 ), YT( 2 )
213 * ..
214 * .. External Subroutines ..
215  EXTERNAL xerbla
216 * ..
217 * .. Intrinsic Functions ..
218  INTRINSIC conjg
219 * ..
220 * .. Executable Statements ..
221 *
222 * Set up indices, arrays for ends
223 *
224  IF( lrows ) THEN
225  iinc = lda
226  inext = 1
227  ELSE
228  iinc = 1
229  inext = lda
230  END IF
231 *
232  IF( lleft ) THEN
233  nt = 1
234  ix = 1 + iinc
235  iy = 2 + lda
236  xt( 1 ) = a( 1 )
237  yt( 1 ) = xleft
238  ELSE
239  nt = 0
240  ix = 1
241  iy = 1 + inext
242  END IF
243 *
244  IF( lright ) THEN
245  iyt = 1 + inext + ( nl-1 )*iinc
246  nt = nt + 1
247  xt( nt ) = xright
248  yt( nt ) = a( iyt )
249  END IF
250 *
251 * Check for errors
252 *
253  IF( nl.LT.nt ) THEN
254  CALL xerbla( 'CLAROT', 4 )
255  RETURN
256  END IF
257  IF( lda.LE.0 .OR. ( .NOT.lrows .AND. lda.LT.nl-nt ) ) THEN
258  CALL xerbla( 'CLAROT', 8 )
259  RETURN
260  END IF
261 *
262 * Rotate
263 *
264 * CROT( NL-NT, A(IX),IINC, A(IY),IINC, C, S ) with complex C, S
265 *
266  DO 10 j = 0, nl - nt - 1
267  tempx = c*a( ix+j*iinc ) + s*a( iy+j*iinc )
268  a( iy+j*iinc ) = -conjg( s )*a( ix+j*iinc ) +
269  $ conjg( c )*a( iy+j*iinc )
270  a( ix+j*iinc ) = tempx
271  10 CONTINUE
272 *
273 * CROT( NT, XT,1, YT,1, C, S ) with complex C, S
274 *
275  DO 20 j = 1, nt
276  tempx = c*xt( j ) + s*yt( j )
277  yt( j ) = -conjg( s )*xt( j ) + conjg( c )*yt( j )
278  xt( j ) = tempx
279  20 CONTINUE
280 *
281 * Stuff values back into XLEFT, XRIGHT, etc.
282 *
283  IF( lleft ) THEN
284  a( 1 ) = xt( 1 )
285  xleft = yt( 1 )
286  END IF
287 *
288  IF( lright ) THEN
289  xright = xt( nt )
290  a( iyt ) = yt( nt )
291  END IF
292 *
293  RETURN
294 *
295 * End of CLAROT
296 *
297  END
clarot
subroutine clarot(LROWS, LLEFT, LRIGHT, NL, C, S, A, LDA, XLEFT, XRIGHT)
Definition: clarot.f:3