LAPACK 3.3.1
Linear Algebra PACKage

clarot.f

Go to the documentation of this file.
00001       SUBROUTINE CLAROT( LROWS, LLEFT, LRIGHT, NL, C, S, A, LDA, XLEFT,
00002      $                   XRIGHT )
00003 *
00004 *  -- LAPACK auxiliary test routine (version 3.1) --
00005 *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
00006 *     November 2006
00007 *
00008 *     .. Scalar Arguments ..
00009       LOGICAL            LLEFT, LRIGHT, LROWS
00010       INTEGER            LDA, NL
00011       COMPLEX            C, S, XLEFT, XRIGHT
00012 *     ..
00013 *     .. Array Arguments ..
00014       COMPLEX            A( * )
00015 *     ..
00016 *
00017 *  Purpose
00018 *  =======
00019 *
00020 *     CLAROT applies a (Givens) rotation to two adjacent rows or
00021 *     columns, where one element of the first and/or last column/row
00022 *     for use on matrices stored in some format other than GE, so
00023 *     that elements of the matrix may be used or modified for which
00024 *     no array element is provided.
00025 *
00026 *     One example is a symmetric matrix in SB format (bandwidth=4), for
00027 *     which UPLO='L':  Two adjacent rows will have the format:
00028 *
00029 *     row j:     *  *  *  *  *  .  .  .  .
00030 *     row j+1:      *  *  *  *  *  .  .  .  .
00031 *
00032 *     '*' indicates elements for which storage is provided,
00033 *     '.' indicates elements for which no storage is provided, but
00034 *     are not necessarily zero; their values are determined by
00035 *     symmetry.  ' ' indicates elements which are necessarily zero,
00036 *      and have no storage provided.
00037 *
00038 *     Those columns which have two '*'s can be handled by SROT.
00039 *     Those columns which have no '*'s can be ignored, since as long
00040 *     as the Givens rotations are carefully applied to preserve
00041 *     symmetry, their values are determined.
00042 *     Those columns which have one '*' have to be handled separately,
00043 *     by using separate variables "p" and "q":
00044 *
00045 *     row j:     *  *  *  *  *  p  .  .  .
00046 *     row j+1:   q  *  *  *  *  *  .  .  .  .
00047 *
00048 *     The element p would have to be set correctly, then that column
00049 *     is rotated, setting p to its new value.  The next call to
00050 *     CLAROT would rotate columns j and j+1, using p, and restore
00051 *     symmetry.  The element q would start out being zero, and be
00052 *     made non-zero by the rotation.  Later, rotations would presumably
00053 *     be chosen to zero q out.
00054 *
00055 *     Typical Calling Sequences: rotating the i-th and (i+1)-st rows.
00056 *     ------- ------- ---------
00057 *
00058 *       General dense matrix:
00059 *
00060 *               CALL CLAROT(.TRUE.,.FALSE.,.FALSE., N, C,S,
00061 *                       A(i,1),LDA, DUMMY, DUMMY)
00062 *
00063 *       General banded matrix in GB format:
00064 *
00065 *               j = MAX(1, i-KL )
00066 *               NL = MIN( N, i+KU+1 ) + 1-j
00067 *               CALL CLAROT( .TRUE., i-KL.GE.1, i+KU.LT.N, NL, C,S,
00068 *                       A(KU+i+1-j,j),LDA-1, XLEFT, XRIGHT )
00069 *
00070 *               [ note that i+1-j is just MIN(i,KL+1) ]
00071 *
00072 *       Symmetric banded matrix in SY format, bandwidth K,
00073 *       lower triangle only:
00074 *
00075 *               j = MAX(1, i-K )
00076 *               NL = MIN( K+1, i ) + 1
00077 *               CALL CLAROT( .TRUE., i-K.GE.1, .TRUE., NL, C,S,
00078 *                       A(i,j), LDA, XLEFT, XRIGHT )
00079 *
00080 *       Same, but upper triangle only:
00081 *
00082 *               NL = MIN( K+1, N-i ) + 1
00083 *               CALL CLAROT( .TRUE., .TRUE., i+K.LT.N, NL, C,S,
00084 *                       A(i,i), LDA, XLEFT, XRIGHT )
00085 *
00086 *       Symmetric banded matrix in SB format, bandwidth K,
00087 *       lower triangle only:
00088 *
00089 *               [ same as for SY, except:]
00090 *                   . . . .
00091 *                       A(i+1-j,j), LDA-1, XLEFT, XRIGHT )
00092 *
00093 *               [ note that i+1-j is just MIN(i,K+1) ]
00094 *
00095 *       Same, but upper triangle only:
00096 *                   . . .
00097 *                       A(K+1,i), LDA-1, XLEFT, XRIGHT )
00098 *
00099 *       Rotating columns is just the transpose of rotating rows, except
00100 *       for GB and SB: (rotating columns i and i+1)
00101 *
00102 *       GB:
00103 *               j = MAX(1, i-KU )
00104 *               NL = MIN( N, i+KL+1 ) + 1-j
00105 *               CALL CLAROT( .TRUE., i-KU.GE.1, i+KL.LT.N, NL, C,S,
00106 *                       A(KU+j+1-i,i),LDA-1, XTOP, XBOTTM )
00107 *
00108 *               [note that KU+j+1-i is just MAX(1,KU+2-i)]
00109 *
00110 *       SB: (upper triangle)
00111 *
00112 *                    . . . . . .
00113 *                       A(K+j+1-i,i),LDA-1, XTOP, XBOTTM )
00114 *
00115 *       SB: (lower triangle)
00116 *
00117 *                    . . . . . .
00118 *                       A(1,i),LDA-1, XTOP, XBOTTM )
00119 *
00120 *  Arguments
00121 *  =========
00122 *
00123 *  LROWS  - LOGICAL
00124 *           If .TRUE., then CLAROT will rotate two rows.  If .FALSE.,
00125 *           then it will rotate two columns.
00126 *           Not modified.
00127 *
00128 *  LLEFT  - LOGICAL
00129 *           If .TRUE., then XLEFT will be used instead of the
00130 *           corresponding element of A for the first element in the
00131 *           second row (if LROWS=.FALSE.) or column (if LROWS=.TRUE.)
00132 *           If .FALSE., then the corresponding element of A will be
00133 *           used.
00134 *           Not modified.
00135 *
00136 *  LRIGHT - LOGICAL
00137 *           If .TRUE., then XRIGHT will be used instead of the
00138 *           corresponding element of A for the last element in the
00139 *           first row (if LROWS=.FALSE.) or column (if LROWS=.TRUE.) If
00140 *           .FALSE., then the corresponding element of A will be used.
00141 *           Not modified.
00142 *
00143 *  NL     - INTEGER
00144 *           The length of the rows (if LROWS=.TRUE.) or columns (if
00145 *           LROWS=.FALSE.) to be rotated.  If XLEFT and/or XRIGHT are
00146 *           used, the columns/rows they are in should be included in
00147 *           NL, e.g., if LLEFT = LRIGHT = .TRUE., then NL must be at
00148 *           least 2.  The number of rows/columns to be rotated
00149 *           exclusive of those involving XLEFT and/or XRIGHT may
00150 *           not be negative, i.e., NL minus how many of LLEFT and
00151 *           LRIGHT are .TRUE. must be at least zero; if not, XERBLA
00152 *           will be called.
00153 *           Not modified.
00154 *
00155 *  C, S   - COMPLEX
00156 *           Specify the Givens rotation to be applied.  If LROWS is
00157 *           true, then the matrix ( c  s )
00158 *                                 ( _  _ )
00159 *                                 (-s  c )  is applied from the left;
00160 *           if false, then the transpose (not conjugated) thereof is
00161 *           applied from the right.  Note that in contrast to the
00162 *           output of CROTG or to most versions of CROT, both C and S
00163 *           are complex.  For a Givens rotation, |C|**2 + |S|**2 should
00164 *           be 1, but this is not checked.
00165 *           Not modified.
00166 *
00167 *  A      - COMPLEX array.
00168 *           The array containing the rows/columns to be rotated.  The
00169 *           first element of A should be the upper left element to
00170 *           be rotated.
00171 *           Read and modified.
00172 *
00173 *  LDA    - INTEGER
00174 *           The "effective" leading dimension of A.  If A contains
00175 *           a matrix stored in GE, HE, or SY format, then this is just
00176 *           the leading dimension of A as dimensioned in the calling
00177 *           routine.  If A contains a matrix stored in band (GB, HB, or
00178 *           SB) format, then this should be *one less* than the leading
00179 *           dimension used in the calling routine.  Thus, if A were
00180 *           dimensioned A(LDA,*) in CLAROT, then A(1,j) would be the
00181 *           j-th element in the first of the two rows to be rotated,
00182 *           and A(2,j) would be the j-th in the second, regardless of
00183 *           how the array may be stored in the calling routine.  [A
00184 *           cannot, however, actually be dimensioned thus, since for
00185 *           band format, the row number may exceed LDA, which is not
00186 *           legal FORTRAN.]
00187 *           If LROWS=.TRUE., then LDA must be at least 1, otherwise
00188 *           it must be at least NL minus the number of .TRUE. values
00189 *           in XLEFT and XRIGHT.
00190 *           Not modified.
00191 *
00192 *  XLEFT  - COMPLEX
00193 *           If LLEFT is .TRUE., then XLEFT will be used and modified
00194 *           instead of A(2,1) (if LROWS=.TRUE.) or A(1,2)
00195 *           (if LROWS=.FALSE.).
00196 *           Read and modified.
00197 *
00198 *  XRIGHT - COMPLEX
00199 *           If LRIGHT is .TRUE., then XRIGHT will be used and modified
00200 *           instead of A(1,NL) (if LROWS=.TRUE.) or A(NL,1)
00201 *           (if LROWS=.FALSE.).
00202 *           Read and modified.
00203 *
00204 *  =====================================================================
00205 *
00206 *     .. Local Scalars ..
00207       INTEGER            IINC, INEXT, IX, IY, IYT, J, NT
00208       COMPLEX            TEMPX
00209 *     ..
00210 *     .. Local Arrays ..
00211       COMPLEX            XT( 2 ), YT( 2 )
00212 *     ..
00213 *     .. External Subroutines ..
00214       EXTERNAL           XERBLA
00215 *     ..
00216 *     .. Intrinsic Functions ..
00217       INTRINSIC          CONJG
00218 *     ..
00219 *     .. Executable Statements ..
00220 *
00221 *     Set up indices, arrays for ends
00222 *
00223       IF( LROWS ) THEN
00224          IINC = LDA
00225          INEXT = 1
00226       ELSE
00227          IINC = 1
00228          INEXT = LDA
00229       END IF
00230 *
00231       IF( LLEFT ) THEN
00232          NT = 1
00233          IX = 1 + IINC
00234          IY = 2 + LDA
00235          XT( 1 ) = A( 1 )
00236          YT( 1 ) = XLEFT
00237       ELSE
00238          NT = 0
00239          IX = 1
00240          IY = 1 + INEXT
00241       END IF
00242 *
00243       IF( LRIGHT ) THEN
00244          IYT = 1 + INEXT + ( NL-1 )*IINC
00245          NT = NT + 1
00246          XT( NT ) = XRIGHT
00247          YT( NT ) = A( IYT )
00248       END IF
00249 *
00250 *     Check for errors
00251 *
00252       IF( NL.LT.NT ) THEN
00253          CALL XERBLA( 'CLAROT', 4 )
00254          RETURN
00255       END IF
00256       IF( LDA.LE.0 .OR. ( .NOT.LROWS .AND. LDA.LT.NL-NT ) ) THEN
00257          CALL XERBLA( 'CLAROT', 8 )
00258          RETURN
00259       END IF
00260 *
00261 *     Rotate
00262 *
00263 *     CROT( NL-NT, A(IX),IINC, A(IY),IINC, C, S ) with complex C, S
00264 *
00265       DO 10 J = 0, NL - NT - 1
00266          TEMPX = C*A( IX+J*IINC ) + S*A( IY+J*IINC )
00267          A( IY+J*IINC ) = -CONJG( S )*A( IX+J*IINC ) +
00268      $                    CONJG( C )*A( IY+J*IINC )
00269          A( IX+J*IINC ) = TEMPX
00270    10 CONTINUE
00271 *
00272 *     CROT( NT, XT,1, YT,1, C, S ) with complex C, S
00273 *
00274       DO 20 J = 1, NT
00275          TEMPX = C*XT( J ) + S*YT( J )
00276          YT( J ) = -CONJG( S )*XT( J ) + CONJG( C )*YT( J )
00277          XT( J ) = TEMPX
00278    20 CONTINUE
00279 *
00280 *     Stuff values back into XLEFT, XRIGHT, etc.
00281 *
00282       IF( LLEFT ) THEN
00283          A( 1 ) = XT( 1 )
00284          XLEFT = YT( 1 )
00285       END IF
00286 *
00287       IF( LRIGHT ) THEN
00288          XRIGHT = XT( NT )
00289          A( IYT ) = YT( NT )
00290       END IF
00291 *
00292       RETURN
00293 *
00294 *     End of CLAROT
00295 *
00296       END
 All Files Functions