Anasazi Version of the Day
Tsqr_CombineFortran.f90
00001 ! @HEADER
00002 ! ***********************************************************************
00003 !
00004 !                 Anasazi: Block Eigensolvers Package
00005 !                 Copyright (2010) Sandia Corporation
00006 !
00007 ! Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 ! license for use of this work by or on behalf of the U.S. Government.
00009 !
00010 ! This library is free software; you can redistribute it and/or modify
00011 ! it under the terms of the GNU Lesser General Public License as
00012 ! published by the Free Software Foundation; either version 2.1 of the
00013 ! License, or (at your option) any later version.
00014 !
00015 ! This library is distributed in the hope that it will be useful, but
00016 ! WITHOUT ANY WARRANTY; without even the implied warranty of
00017 ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 ! Lesser General Public License for more details.
00019 !
00020 ! You should have received a copy of the GNU Lesser General Public
00021 ! License along with this library; if not, write to the Free Software
00022 ! Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 ! USA
00024 ! Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00025 !
00026 ! ***********************************************************************
00027 ! @HEADER
00028 
00029 module TsqrCombine
00030   use TsqrHouseholderReflector, only : DLARFP_wrapper, SLARFP_wrapper, ZLARFP_wrapper, CLARFP_wrapper
00031 
00032   implicit none
00033 
00034   interface 
00035      subroutine DGEMV(TRANS,M,N,ALPHA,A,LDA,X,INCX,BETA,Y,INCY)
00036        real (8), intent(in)    :: ALPHA,BETA
00037        integer, intent(in)     :: INCX,INCY,LDA,M,N
00038        character, intent(in)   :: TRANS
00039        real (8), intent(in)    :: A(LDA,*), X(*)
00040        real (8), intent(inout) :: Y(*)
00041      end subroutine DGEMV
00042 
00043      subroutine ZGEMV(TRANS,M,N,ALPHA,A,LDA,X,INCX,BETA,Y,INCY)
00044        complex (8), intent(in)    :: ALPHA,BETA
00045        integer, intent(in)        :: INCX,INCY,LDA,M,N
00046        character, intent(in)      :: TRANS
00047        complex (8), intent(in)    :: A(LDA,*), X(*)
00048        complex (8), intent(inout) :: Y(*)
00049      end subroutine ZGEMV
00050 
00051      subroutine SGEMV(TRANS,M,N,ALPHA,A,LDA,X,INCX,BETA,Y,INCY)
00052        real (4), intent(in)    :: ALPHA,BETA
00053        integer, intent(in)     :: INCX,INCY,LDA,M,N
00054        character, intent(in)   :: TRANS
00055        real (4), intent(in)    :: A(LDA,*), X(*)
00056        real (4), intent(inout) :: Y(*)
00057      end subroutine SGEMV
00058 
00059      subroutine CGEMV(TRANS,M,N,ALPHA,A,LDA,X,INCX,BETA,Y,INCY)
00060        complex (4), intent(in)    :: ALPHA,BETA
00061        integer, intent(in)        :: INCX,INCY,LDA,M,N
00062        character, intent(in)      :: TRANS
00063        complex (4), intent(in)    :: A(LDA,*), X(*)
00064        complex (4), intent(inout) :: Y(*)
00065      end subroutine CGEMV
00066      
00067      subroutine DGER(M,N,ALPHA,X,INCX,Y,INCY,A,LDA)
00068        real (8), intent(in)    :: ALPHA
00069        integer, intent(in)     :: INCX,INCY,LDA,M,N
00070        real (8), intent(inout) :: A(LDA,*)
00071        real (8), intent(in)    :: X(*),Y(*)
00072      end subroutine DGER
00073 
00074      subroutine ZGERC(M,N,ALPHA,X,INCX,Y,INCY,A,LDA)
00075        complex (8), intent(in)    :: ALPHA
00076        integer, intent(in)        :: INCX,INCY,LDA,M,N
00077        complex (8), intent(inout) :: A(LDA,*)
00078        complex (8), intent(in)    :: X(*),Y(*)
00079      end subroutine ZGERC
00080 
00081      subroutine SGER(M,N,ALPHA,X,INCX,Y,INCY,A,LDA)
00082        real (4), intent(in)    :: ALPHA
00083        integer, intent(in)     :: INCX,INCY,LDA,M,N
00084        real (4), intent(inout) :: A(LDA,*)
00085        real (4), intent(in)    :: X(*),Y(*)
00086      end subroutine SGER
00087 
00088      subroutine CGERC(M,N,ALPHA,X,INCX,Y,INCY,A,LDA)
00089        complex (4), intent(in)    :: ALPHA
00090        integer, intent(in)        :: INCX,INCY,LDA,M,N
00091        complex (4), intent(inout) :: A(LDA,*)
00092        complex (4), intent(in)    :: X(*),Y(*)
00093      end subroutine CGERC
00094   end interface
00095 
00096 contains
00097 
00098   ! Apply the Q factor stored in [R; A] to [C_top; C_bot].  The C
00099   ! blocks are allowed, but not required, to have different leading
00100   ! dimensions (ldc_top resp. ldc_bottom).  R is upper triangular, so
00101   ! we do not need it; the Householder reflectors representing the Q
00102   ! factor are stored compactly in A (specifically, in all of A, not
00103   ! just the lower triangle).
00104   !
00105   ! In the "sequential under parallel" version of TSQR, this function
00106   ! belongs to the sequential part (i.e., operating on cache blocks on
00107   ! a single processor).
00108   !
00109   ! This routine used to be called "tsqr_apply_inner2" (the original
00110   ! tsqr_apply_inner assumed that C_top and C_bot had the same leading
00111   ! dimension, but was otherwise the same).  After that, it was called
00112   ! "tsqr_apply_inner".  Now we have different function names, based
00113   ! on input scalar datatypes: s_apply_inner, d_apply_inner, etc.  We
00114   ! also removed the ISO_C_BINDING dependency, since that doesn't work
00115   ! with older compilers.
00116   !
00117   ! trans [in]     'N' means apply Q, anything else (such as 'T') 
00118   !                means apply Q^T
00119   ! m [in]         number of rows of A
00120   ! ncols_C [in]   number of columns of [C_top; C_bot]
00121   ! ncols_Q [in]   number of columns of [R; A]
00122   ! A [in]         m by ncols_Q matrix, in which the Householder 
00123   !                reflectors representing the Q factor are stored
00124   ! lda [in]       leading dimension of A
00125   ! tau [in]       array of length ncols_Q, storing the scaling factors 
00126   !                for the Householder reflectors representing Q 
00127   ! C_top [inout]  ncols_Q by ncols_C matrix
00128   ! ldc_top [in]   leading dimension of C_top
00129   ! C_bot [inout]  m by ncols_C matrix
00130   ! ldc_bot [in]   leading dimension of C_bot
00131   ! work [out]     workspace array of length ncols_C
00132   ! 
00133   subroutine d_apply_inner( trans, m, ncols_C, ncols_Q, A, lda, tau, &
00134        C_top, ldc_top, C_bot, ldc_bot, work )
00135     implicit none
00136     
00137     character, intent(in)        :: trans
00138     integer, intent(in), value   :: m, ncols_Q, ncols_C, lda, ldc_top, ldc_bot
00139     real(8), intent(in)          :: A(lda,ncols_Q)
00140     real(8), intent(inout)       :: C_top(ldc_top,ncols_C), C_bot(ldc_bot,ncols_C)
00141     real(8), intent(in)          :: tau(ncols_Q)
00142     real(8), intent(out), target :: work(ncols_C)
00143 
00144     real (8), pointer :: y(:)
00145     integer           :: j, k
00146 
00147     y => work(1:ncols_C)
00148     y = 0
00149     if (trans == 'N' .or. trans == 'n') then
00150        do k = ncols_Q, 1, -1
00151           ! The m bottom elements of the length m+1 "sparse" Householder
00152           ! reflector are stored in A(1:m,k).  The topmost element is
00153           ! implicit and is 1.0.  If we want to exploit its sparsity, we
00154           ! can't use DLARF to apply it to [C_top; C_bot]; we have to
00155           ! roll our own routine, which you see here.
00156           
00157           ! $y^T := A(1:m,k)^T C_bot(1:m, 1:ncols_C)$, so $y := C_bot(1:m, 1:ncols_C)^T A(1:m, k)$.
00158           !call DGEMV( 'T', m, ncols_C, 1.0d0, C_bot(1:m, 1:ncols_C), ldc, A(1:m, k), 1, 0.0d0, y(1:ncols_C), 1 )
00159           call DGEMV( 'T', m, ncols_C, 1.0d0, C_bot, ldc_bot, A(1, k), 1, 0.0d0, y, 1 )
00160           
00161           ! $y^T := y^T + C_top(k, 1:ncols_C)$
00162           do j = 1, ncols_C
00163              y(j) = y(j) + C_top(k, j)
00164           end do
00165           
00166           ! Update C_top(k, 1:ncols_C)
00167           do j = 1, ncols_C
00168              C_top(k, j) = C_top(k, j) - tau(k) * y(j)
00169           end do
00170           
00171           ! Update C_bot(1:m, 1:ncols_C)
00172           !call DGER( m, ncols_C, -tau(k), A(1:m,k), 1, y(1:ncols_C), 1, C_bot(1:m, 1:ncols_C), ldc )
00173           call DGER( m, ncols_C, -tau(k), A(1,k), 1, y, 1, C_bot, ldc_bot )
00174        end do
00175     else
00176        do k = 1, ncols_Q
00177           ! The m bottom elements of the length m+1 "sparse" Householder
00178           ! reflector are stored in A(1:m,k).  The topmost element is
00179           ! implicit and is 1.0.  If we want to exploit its sparsity, we
00180           ! can't use DLARF to apply it to [C_top; C_bot]; we have to
00181           ! roll our own routine, which you see here.
00182           
00183           ! $y^T := A(1:m,k)^T C_bot(1:m, 1:ncols_C)$, so $y := C_bot(1:m, 1:ncols_C)^T A(1:m, k)$.
00184           !call DGEMV( 'T', m, ncols_C, 1.0d0, C_bot(1:m, 1:ncols_C), ldc_bot, A(1:m, k), 1, 0.0d0, y(1:ncols_C), 1 )
00185           call DGEMV( 'T', m, ncols_C, 1.0d0, C_bot, ldc_bot, A(1, k), 1, 0.0d0, y, 1 )
00186           
00187           ! $y^T := y^T + C_top(k, 1:ncols_C)$
00188           do j = 1, ncols_C
00189              y(j) = y(j) + C_top(k, j)
00190           end do
00191           
00192           ! Update C_top(k, 1:ncols_C)
00193           do j = 1, ncols_C
00194              C_top(k, j) = C_top(k, j) - tau(k) * y(j)
00195           end do
00196           
00197           ! Update C_bot(1:m, 1:ncols_C)
00198           !call DGER( m, ncols_C, -tau(k), A(1:m,k), 1, y(1:ncols_C), 1, C_bot(1:m, 1:ncols_C), ldc_bot )
00199           call DGER( m, ncols_C, -tau(k), A(1,k), 1, y, 1, C_bot(1, 1), ldc_bot )
00200        end do
00201     end if
00202   end subroutine d_apply_inner
00203 
00204 
00205   subroutine z_apply_inner( trans, m, ncols_C, ncols_Q, A, lda, tau, &
00206        C_top, ldc_top, C_bot, ldc_bot, work )
00207     implicit none
00208     
00209     character, intent(in)           :: trans
00210     integer, intent(in), value      :: m, ncols_Q, ncols_C, lda, ldc_top, ldc_bot
00211     complex(8), intent(in)          :: A(lda,ncols_Q)
00212     complex(8), intent(inout)       :: C_top(ldc_top,ncols_C), C_bot(ldc_bot,ncols_C)
00213     complex(8), intent(in)          :: tau(ncols_Q)
00214     complex(8), intent(out)         :: work(ncols_C)
00215 
00216     complex(8)                      :: ZERO, ONE, tau_k
00217     integer                         :: j, k, k_first, k_second, k_step
00218     logical                         :: no_trans
00219 
00220     ZERO = complex( 0.0d0, 0.0d0 )
00221     ONE = complex( 1.0d0, 0.0d0 )
00222     do k = 1, ncols_C
00223        work(k) = ZERO
00224     end do
00225 
00226     no_trans = (trans == 'N' .or. trans == 'n')
00227     if (no_trans) then
00228        k_first = ncols_Q
00229        k_second = 1
00230        k_step = -1
00231     else 
00232        k_first = 1
00233        k_second = ncols_Q
00234        k_step = 1
00235     end if
00236 
00237     do k = k_first, k_second, k_step
00238        if (no_trans) then
00239           tau_k = tau(k)
00240        else
00241           tau_k = conjg( tau(k) )
00242        end if
00243 
00244        call ZGEMV( 'Conjugate transpose', m, ncols_C, ONE, C_bot, ldc_bot, A(1, k), 1, ZERO, work, 1 )
00245        do j = 1, ncols_C
00246           work(j) = work(j) + conjg( C_top(k, j) )
00247        end do
00248        do j = 1, ncols_C
00249           C_top(k, j) = C_top(k, j) - tau_k * work(j)
00250        end do
00251        call ZGERC( m, ncols_C, -tau_k, A(1,k), 1, work, 1, C_bot(1, 1), ldc_bot )
00252     end do
00253   end subroutine z_apply_inner
00254 
00255 
00256 
00257   subroutine s_apply_inner( trans, m, ncols_C, ncols_Q, A, lda, tau, &
00258        C_top, ldc_top, C_bot, ldc_bot, work )
00259     implicit none
00260     
00261     character, intent(in)        :: trans
00262     integer, intent(in), value   :: m, ncols_Q, ncols_C, lda, ldc_top, ldc_bot
00263     real(4), intent(in)          :: A(lda,ncols_Q)
00264     real(4), intent(inout)       :: C_top(ldc_top,ncols_C), C_bot(ldc_bot,ncols_C)
00265     real(4), intent(in)          :: tau(ncols_Q)
00266     real(4), intent(out), target :: work(ncols_C)
00267 
00268     real(4), pointer :: y(:)
00269     integer          :: j, k
00270 
00271     y => work(1:ncols_C)
00272     y = 0
00273     if (trans == 'N' .or. trans == 'n') then
00274        do k = ncols_Q, 1, -1
00275           ! The m bottom elements of the length m+1 "sparse"
00276           ! Householder reflector are stored in A(1:m,k).  The topmost
00277           ! element is implicit and is 1.0.  If we want to exploit its
00278           ! sparsity, we can't use SLARF to apply it to [C_top;
00279           ! C_bot]; we have to roll our own routine, which you see
00280           ! here.
00281           !
00282           ! $y^T := A(1:m,k)^T C_bot(1:m, 1:ncols_C)$, so $y := C_bot(1:m, 1:ncols_C)^T A(1:m, k)$.
00283           call SGEMV( 'T', m, ncols_C, 1.0, C_bot, ldc_bot, A(1, k), 1, 0.0, y, 1 )
00284           
00285           ! $y^T := y^T + C_top(k, 1:ncols_C)$
00286           do j = 1, ncols_C
00287              y(j) = y(j) + C_top(k, j)
00288           end do
00289           
00290           ! Update C_top(k, 1:ncols_C)
00291           do j = 1, ncols_C
00292              C_top(k, j) = C_top(k, j) - tau(k) * y(j)
00293           end do
00294           
00295           ! Update C_bot(1:m, 1:ncols_C)
00296           call SGER( m, ncols_C, -tau(k), A(1,k), 1, y, 1, C_bot, ldc_bot )
00297        end do
00298     else
00299        do k = 1, ncols_Q
00300           ! The m bottom elements of the length m+1 "sparse"
00301           ! Householder reflector are stored in A(1:m,k).  The topmost
00302           ! element is implicit and is 1.0.  If we want to exploit its
00303           ! sparsity, we can't use SLARF to apply it to [C_top;
00304           ! C_bot]; we have to roll our own routine, which you see
00305           ! here.
00306           !
00307           ! $y^T := A(1:m,k)^T C_bot(1:m, 1:ncols_C)$, so $y := C_bot(1:m, 1:ncols_C)^T A(1:m, k)$.
00308           call SGEMV( 'T', m, ncols_C, 1.0, C_bot, ldc_bot, A(1, k), 1, 0.0, y, 1 )
00309           
00310           ! $y^T := y^T + C_top(k, 1:ncols_C)$
00311           do j = 1, ncols_C
00312              y(j) = y(j) + C_top(k, j)
00313           end do
00314           
00315           ! Update C_top(k, 1:ncols_C)
00316           do j = 1, ncols_C
00317              C_top(k, j) = C_top(k, j) - tau(k) * y(j)
00318           end do
00319           
00320           ! Update C_bot(1:m, 1:ncols_C)
00321           call SGER( m, ncols_C, -tau(k), A(1,k), 1, y, 1, C_bot(1, 1), ldc_bot )
00322        end do
00323     end if
00324   end subroutine s_apply_inner
00325 
00326 
00327 
00328   subroutine c_apply_inner( trans, m, ncols_C, ncols_Q, A, lda, tau, &
00329        C_top, ldc_top, C_bot, ldc_bot, work )
00330     implicit none
00331     
00332     character, intent(in)      :: trans
00333     integer, intent(in), value :: m, ncols_Q, ncols_C, lda, ldc_top, ldc_bot
00334     complex(4), intent(in)     :: A(lda,ncols_Q)
00335     complex(4), intent(inout)  :: C_top(ldc_top,ncols_C), C_bot(ldc_bot,ncols_C)
00336     complex(4), intent(in)     :: tau(ncols_Q)
00337     complex(4), intent(out)    :: work(ncols_C)
00338 
00339     complex(4)                 :: ZERO, ONE, tau_k
00340     integer                    :: j, k, k_first, k_second, k_step
00341     logical                    :: no_trans
00342 
00343     ZERO = complex( 0.0e0, 0.0e0 )
00344     ONE = complex( 1.0e0, 0.0e0 )
00345     do k = 1, ncols_C
00346        work(k) = ZERO
00347     end do
00348 
00349     no_trans = (trans == 'N' .or. trans == 'n')
00350     if (no_trans) then
00351        k_first = ncols_Q
00352        k_second = 1
00353        k_step = -1
00354     else 
00355        k_first = 1
00356        k_second = ncols_Q
00357        k_step = 1
00358     end if
00359 
00360     do k = k_first, k_second, k_step
00361        if (no_trans) then
00362           tau_k = tau(k)
00363        else
00364           tau_k = conjg( tau(k) )
00365        end if
00366 
00367        call CGEMV( 'Conjugate transpose', m, ncols_C, ONE, C_bot, ldc_bot, &
00368             A(1, k), 1, ZERO, work, 1 )
00369        do j = 1, ncols_C
00370           work(j) = work(j) + conjg( C_top(k, j) )
00371        end do
00372        do j = 1, ncols_C
00373           C_top(k, j) = C_top(k, j) - tau_k * work(j)
00374        end do
00375        call CGERC( m, ncols_C, -tau_k, A(1,k), 1, work, 1, C_bot(1, 1), ldc_bot )
00376     end do
00377   end subroutine c_apply_inner
00378 
00379 
00380   ! Perform one "inner" QR factorization step of sequential / parallel
00381   ! TSQR.  (In either case, only one processor calls this routine.)
00382   !
00383   ! In the "sequential under parallel" version of TSQR, this function
00384   ! belongs to the sequential part (i.e., operating on cache blocks on
00385   ! a single processor).  Only the first cache block $A_0$ is factored
00386   ! as $Q_0 R_0 = A_0$ (see tsqr_factor_first); subsequent cache blocks
00387   ! $A_k$ are factored using this routine, which combines them with 
00388   ! $R_{k-1}$.
00389   !
00390   ! Here is the matrix to factor:
00391   ! \[
00392   ! \begin{pmatrix}
00393   ! R_{k-1} \\      % $A_{k-1}$ is $m_{k-1} \times n$ with $m_{k-1} \geq n$
00394   ! A_k     \\      % $m_k \times n$ with $m_k \geq n$
00395   ! \end{pmatrix}
00396   ! \]
00397   !
00398   ! Since $R_{k-1}$ is n by n upper triangular, we can store the
00399   ! Householder reflectors (representing the Q factor of $[R_{k-1};
00400   ! A_k]$) entirely in $A_k$ (specifically, in all of $A_k$, not just
00401   ! below the diagonal).
00402   !
00403   ! m [in]     Number of rows in the "bottom" block to factor.  The number of
00404   !            rows in the top block doesn't matter, given the assumptions 
00405   !            above, as long as $m_{k-1} \geq n$.
00406   ! n [in]     Number of columns (same in both blocks)
00407   ! R [inout]  "Top" upper triangular n by n block $R_{k-1}$.
00408   !            Overwritten with the new R factor $R_k$ of $[R_{k-1}; A_k]$.
00409   ! ldr [in]   Leading dimension of R
00410   ! A [inout]  "Bottom" dense m by n block $A_k$.  Overwritten with the 
00411   !            Householder reflectors representing the Q factor of 
00412   !            $[R_{k-1}; A_k]$.
00413   ! tau [out]  Scaling factors of the Householder reflectors.  Corresponds 
00414   !            to the output of DGEQRF of the same name.
00415   ! work [out] Workspace (length >= n; don't need lwork or workspace query)
00416   !
00417   subroutine d_factor_inner( m, n, R, ldr, A, lda, tau, work )
00418     implicit none
00419 
00420     integer, value, intent(in)    :: m, n, ldr, lda
00421     real(8), intent(inout)        :: R(ldr,n)
00422     real(8), intent(inout)        :: A(lda,n)
00423     real(8), intent(out)          :: tau(n)
00424     ! "target" means it's legitimate to have pointers alias this
00425     ! array, internal to this function.
00426     real(8), intent(out), target  :: work(n)
00427 
00428     real(8), pointer  :: y(:)
00429     integer           :: j, k
00430 
00431     y => work(1:n)
00432     y = 0
00433     do k = 1, n-1
00434        ! Form the "sparse" Householder reflector, so that the diagonal
00435        ! elements of the R factor are nonnegative.
00436        ! call DLARFP( m + 1, R(k,k), A(1:m,k), 1, tau(k) )
00437        call DLARFP_wrapper( m + 1, R(k,k), A(1,k), 1, tau(k) )
00438 
00439        ! $y^T := A(1:m,k)^T A(1:m, k+1:n)$, so $y := A(1:m, k+1:n)^T A(1:m, k)$.
00440        ! BEGIN mfh 07 June 2009
00441        ! Another segfault here, which could mean that the use of subarrays 
00442        ! must somehow involve allocating and copying data -- not what we want.
00443        !call DGEMV( 'T', m, n-k, 1.0d0, A(1:m, k+1:n), lda, A(1:m, k), 1, 0.0d0, y(1:n-k), 1 )
00444        call DGEMV( 'T', m, n-k, 1.0d0, A(1, k+1), lda, A(1, k), 1, 0.0d0, y, 1 )
00445 
00446        ! $y^T := y^T + R(k, k+1:n)$
00447        do j = k+1, n
00448           y(j-k) = y(j-k) + R(k, j)
00449        end do
00450 
00451        ! Update R(k, k+1:n)
00452        do j = k+1, n
00453           R(k, j) = R(k, j) - tau(k) * y(j-k)
00454        end do
00455 
00456        ! Update A(1:m, k+1:n)
00457        ! BEGIN mfh 07 June 2009
00458        ! free() was triggering a segfault here, which means the use of subarrays 
00459        ! must somehow involve allocating and copying data -- not what we want.
00460        !call DGER( m, n-k, -tau(k), A(1:m,k), 1, y(1:n-k), 1, A(1:m, k+1:n), lda )
00461        call DGER( m, n-k, -tau(k), A(1,k), 1, y, 1, A(1, k+1), lda )
00462        ! END mfh 07 June 2009
00463     end do
00464 
00465     ! Compute the Householder reflector for the last column.  This
00466     ! last iteration doesn't require an update of the trailing matrix,
00467     ! because there is no trailing matrix left!
00468     !call DLARFP( m + 1, R(n,n), A(1:m,n), 1, tau(n) )
00469     call DLARFP_wrapper( m + 1, R(n,n), A(1,n), 1, tau(n) )
00470   end subroutine d_factor_inner
00471 
00472 
00473   subroutine z_factor_inner( m, n, R, ldr, A, lda, tau, work )
00474     implicit none
00475 
00476     integer, value, intent(in) :: m, n, ldr, lda
00477     complex (8), intent(inout) :: R(ldr,n)
00478     complex (8), intent(inout) :: A(lda,n)
00479     complex (8), intent(out)   :: tau(n)
00480     complex (8), intent(out)   :: work(n)
00481 
00482     complex (8)                :: ZERO, ONE
00483     integer                    :: j, k
00484 
00485     ZERO = complex( 0.0d0, 0.0d0 )
00486     ONE = complex( 1.0d0, 0.0d0 )
00487     do k = 1, n
00488        work(k) = ZERO
00489     end do
00490 
00491     do k = 1, n-1
00492        ! Form the "sparse" Householder reflector, so that the diagonal
00493        ! elements of the R factor are nonnegative.
00494        call ZLARFP_wrapper( m + 1, R(k,k), A(1,k), 1, tau(k) )
00495 
00496        ! $y^* := A(1:m,k)^* A(1:m, k+1:n)$, so $y := A(1:m, k+1:n)^* A(1:m, k)$.
00497        call ZGEMV( 'Conjugate transpose', m, n-k, ONE, A(1, k+1), lda, &
00498             A(1, k), 1, ZERO, work, 1 )
00499 
00500        ! $y^* := y^* + R(k, k+1:n)$, so $y := y + R(k,k+1:n)^*$.
00501        do j = k+1, n
00502           work(j-k) = work(j-k) + conjg( R(k, j) )
00503        end do
00504 
00505        ! Update R(k, k+1:n)
00506        do j = k+1, n
00507           R(k, j) = R(k, j) - tau(k) * work(j-k)
00508        end do
00509 
00510        ! $A(1:m, k+1:n) := A(1:m, k+1:n) + \alpha \cdot x \cdot y^*$
00511        call ZGERC( m, n-k, -tau(k), A(1,k), 1, work, 1, A(1, k+1), lda )
00512     end do
00513 
00514     ! Compute the Householder reflector for the last column.
00515     call ZLARFP_wrapper( m + 1, R(n,n), A(1,n), 1, tau(n) )
00516   end subroutine z_factor_inner
00517 
00518 
00519   subroutine s_factor_inner( m, n, R, ldr, A, lda, tau, work )
00520     implicit none
00521 
00522     integer, value, intent(in)    :: m, n, ldr, lda
00523     real(4), intent(inout)        :: R(ldr,n)
00524     real(4), intent(inout)        :: A(lda,n)
00525     real(4), intent(out)          :: tau(n)
00526     ! "target" means it's legitimate to have pointers alias this
00527     ! array, internal to this function.
00528     real(4), intent(out), target  :: work(n)
00529 
00530     real(4), pointer  :: y(:)
00531     integer           :: j, k
00532 
00533     y => work(1:n)
00534     y = 0
00535     do k = 1, n-1
00536        ! Form the "sparse" Householder reflector, so that the diagonal
00537        ! elements of the R factor are nonnegative.
00538        call SLARFP_WRAPPER( m + 1, R(k,k), A(1,k), 1, tau(k) )
00539 
00540        ! $y^T := A(1:m,k)^T A(1:m, k+1:n)$, so $y := A(1:m, k+1:n)^T A(1:m, k)$.
00541        call SGEMV( 'T', m, n-k, 1.0, A(1, k+1), lda, A(1, k), 1, 0.0, y, 1 )
00542 
00543        ! $y^T := y^T + R(k, k+1:n)$
00544        do j = k+1, n
00545           y(j-k) = y(j-k) + R(k, j)
00546        end do
00547 
00548        ! Update R(k, k+1:n)
00549        do j = k+1, n
00550           R(k, j) = R(k, j) - tau(k) * y(j-k)
00551        end do
00552 
00553        ! Update A(1:m, k+1:n)
00554        call SGER( m, n-k, -tau(k), A(1,k), 1, y, 1, A(1, k+1), lda )
00555     end do
00556 
00557     ! Compute the Householder reflector for the last column.  This
00558     ! last iteration doesn't require an update of the trailing matrix,
00559     ! because there is no trailing matrix left!
00560     call SLARFP_WRAPPER( m + 1, R(n,n), A(1,n), 1, tau(n) )
00561   end subroutine s_factor_inner
00562 
00563 
00564   subroutine c_factor_inner( m, n, R, ldr, A, lda, tau, work )
00565     implicit none
00566 
00567     integer, value, intent(in) :: m, n, ldr, lda
00568     complex (4), intent(inout) :: R(ldr,n)
00569     complex (4), intent(inout) :: A(lda,n)
00570     complex (4), intent(out)   :: tau(n)
00571     complex (4), intent(out)   :: work(n)
00572 
00573     complex (4)                :: ZERO, ONE
00574     integer                    :: j, k
00575 
00576     ZERO = complex( 0.0e0, 0.0e0 )
00577     ONE = complex( 1.0e0, 0.0e0 )
00578     do k = 1, n
00579        work(k) = ZERO
00580     end do
00581 
00582     do k = 1, n-1
00583        ! Form the "sparse" Householder reflector, so that the diagonal
00584        ! elements of the R factor are nonnegative.
00585        call CLARFP_wrapper( m + 1, R(k,k), A(1,k), 1, tau(k) )
00586 
00587        ! $y^* := A(1:m,k)^* A(1:m, k+1:n)$, so $y := A(1:m, k+1:n)^* A(1:m, k)$.
00588        call CGEMV( 'Conjugate transpose', m, n-k, ONE, A(1, k+1), lda, A(1, k), 1, ZERO, work, 1 )
00589 
00590        ! $y^* := y^* + R(k, k+1:n)$, so $y := y + R(k,k+1:n)^*$.
00591        do j = k+1, n
00592           work(j-k) = work(j-k) + conjg( R(k, j) )
00593        end do
00594 
00595        ! Update R(k, k+1:n)
00596        do j = k+1, n
00597           R(k, j) = R(k, j) - tau(k) * work(j-k)
00598        end do
00599 
00600        ! $A(1:m, k+1:n) := A(1:m, k+1:n) + \alpha \cdot x \cdot y^*$
00601        call CGERC( m, n-k, -tau(k), A(1,k), 1, work, 1, A(1, k+1), lda )
00602     end do
00603 
00604     ! Compute the Householder reflector for the last column.
00605     call CLARFP_wrapper( m + 1, R(n,n), A(1,n), 1, tau(n) )
00606   end subroutine c_factor_inner
00607 
00608 
00609 
00610   ! Compute QR factorization of [R_top; R_bot].  Store resulting R
00611   ! factor in R_top and Householder reflectors in R_bot.
00612   !
00613   ! n [in]         Number of rows and columns of each of R_top and R_bot
00614   ! R_top [inout]  n by n upper triangular matrix 
00615   ! ldr_top [in]   Leading dimension of R_top
00616   ! R_bot [inout]  n by n upper triangular matrix 
00617   ! ldr_bot [in]   Leading dimension of R_bot
00618   ! tau [out]      Scaling factors for Householder reflectors
00619   ! work [out]     Workspace array (of length >= n)
00620   !
00621   subroutine d_factor_pair( n, R_top, ldr_top, R_bot, ldr_bot, tau, work )
00622     implicit none
00623 
00624     integer, value, intent(in)    :: n, ldr_top, ldr_bot
00625     real(8), intent(inout)        :: R_top(ldr_top,n), R_bot(ldr_bot,n)
00626     real(8), intent(out)          :: tau(n)
00627     ! "target" means it's legitimate to have pointers alias this
00628     ! array, internal to this function.
00629     real(8), intent(out), target  :: work(n)
00630 
00631     real(8), pointer  :: y(:)
00632     integer           :: j, k
00633 
00634     y => work(1:n)
00635     y = 0
00636     do k = 1, n-1
00637        ! Form the "sparse" Householder reflector, so that the diagonal
00638        ! elements of the R factor are nonnegative.  Length of this
00639        ! reflector is k+1, including the one top element (on the
00640        ! diagonal of R_top) and k bottom elements (above and including
00641        ! the diagonal of R_bot).
00642        call DLARFP_wrapper( k + 1, R_top(k,k), R_bot(1,k), 1, tau(k) )
00643 
00644        ! $y^T := R_bot(1:k,k)^T R_bot(1:k, k+1:n)$, so $y := R_bot(1:k, k+1:n)^T R_bot(1:k, k)$.
00645        call DGEMV( 'T', k, n-k, 1.0d0, R_bot(1, k+1), ldr_bot, R_bot(1, k), 1, 0.0d0, y, 1 )
00646 
00647        ! $y^T := y^T + R_top(k, k+1:n)$
00648        do j = k+1, n
00649           y(j-k) = y(j-k) + R_top(k, j)
00650        end do
00651 
00652        ! Update R_top(k, k+1:n)
00653        do j = k+1, n
00654           R_top(k, j) = R_top(k, j) - tau(k) * y(j-k)
00655        end do
00656 
00657        ! Update R_bot(1:k, k+1:n)
00658        call DGER( k, n-k, -tau(k), R_bot(1,k), 1, y, 1, R_bot(1, k+1), ldr_bot )
00659     end do
00660 
00661     ! Compute the Householder reflector for the last column.  This
00662     ! last iteration doesn't require an update of the trailing matrix,
00663     ! because there is no trailing matrix left!
00664     call DLARFP_wrapper( n + 1, R_top(n,n), R_bot(1,n), 1, tau(n) )
00665   end subroutine d_factor_pair
00666 
00667 
00668   subroutine z_factor_pair( n, R_top, ldr_top, R_bot, ldr_bot, tau, work )
00669     implicit none
00670 
00671     integer, value, intent(in) :: n, ldr_top, ldr_bot
00672     complex (8), intent(inout) :: R_top(ldr_top,n), R_bot(ldr_bot,n)
00673     complex (8), intent(out)   :: tau(n)
00674     complex (8), intent(out)   :: work(n)
00675 
00676     complex (8)                :: ZERO, ONE
00677     integer                    :: j, k
00678 
00679     ZERO = complex( 0.0d0, 0.0d0 )
00680     ONE = complex( 1.0d0, 0.0d0 )
00681     do k = 1, n
00682        work(k) = ZERO
00683     end do
00684 
00685     do k = 1, n-1
00686        call ZLARFP_wrapper( k + 1, R_top(k,k), R_bot(1,k), 1, tau(k) )
00687        call ZGEMV( 'Conjugate transpose', k, n-k, ONE, &
00688             R_bot(1, k+1), ldr_bot, &
00689             R_bot(1, k), 1, ZERO, work, 1 )
00690        do j = k+1, n
00691           work(j-k) = work(j-k) + conjg( R_top(k, j) )
00692        end do
00693        do j = k+1, n
00694           R_top(k, j) = R_top(k, j) - tau(k) * work(j-k)
00695        end do
00696        call ZGERC( k, n-k, -tau(k), R_bot(1,k), 1, work, 1, &
00697             R_bot(1, k+1), ldr_bot )
00698     end do
00699 
00700     call ZLARFP_wrapper( n + 1, R_top(n,n), R_bot(1,n), 1, tau(n) )
00701   end subroutine z_factor_pair
00702 
00703 
00704   subroutine s_factor_pair( n, R_top, ldr_top, R_bot, ldr_bot, tau, work )
00705     implicit none
00706 
00707     integer, value, intent(in)    :: n, ldr_top, ldr_bot
00708     real(4), intent(inout)        :: R_top(ldr_top,n), R_bot(ldr_bot,n)
00709     real(4), intent(out)          :: tau(n)
00710     ! "target" means it's legitimate to have pointers alias this
00711     ! array, internal to this function.
00712     real(4), intent(out), target  :: work(n)
00713 
00714     real(4), pointer  :: y(:)
00715     integer           :: j, k
00716 
00717     y => work(1:n)
00718     y = 0
00719     do k = 1, n-1
00720        ! Form the "sparse" Householder reflector, so that the diagonal
00721        ! elements of the R factor are nonnegative.  Length of this
00722        ! reflector is k+1, including the one top element (on the
00723        ! diagonal of R_top) and k bottom elements (above and including
00724        ! the diagonal of R_bot).
00725        call SLARFP_WRAPPER( k + 1, R_top(k,k), R_bot(1,k), 1, tau(k) )
00726 
00727        ! $y^T := R_bot(1:k,k)^T R_bot(1:k, k+1:n)$, so $y := R_bot(1:k, k+1:n)^T R_bot(1:k, k)$.
00728        call SGEMV( 'T', k, n-k, 1.0, R_bot(1, k+1), ldr_bot, R_bot(1, k), 1, 0.0, y, 1 )
00729 
00730        ! $y^T := y^T + R_top(k, k+1:n)$
00731        do j = k+1, n
00732           y(j-k) = y(j-k) + R_top(k, j)
00733        end do
00734 
00735        ! Update R_top(k, k+1:n)
00736        do j = k+1, n
00737           R_top(k, j) = R_top(k, j) - tau(k) * y(j-k)
00738        end do
00739 
00740        ! Update R_bot(1:k, k+1:n)
00741        call SGER( k, n-k, -tau(k), R_bot(1,k), 1, y, 1, R_bot(1, k+1), ldr_bot )
00742     end do
00743 
00744     ! Compute the Householder reflector for the last column.  This
00745     ! last iteration doesn't require an update of the trailing matrix,
00746     ! because there is no trailing matrix left!
00747     call SLARFP_WRAPPER( n + 1, R_top(n,n), R_bot(1,n), 1, tau(n) )
00748   end subroutine s_factor_pair
00749 
00750 
00751   subroutine c_factor_pair( n, R_top, ldr_top, R_bot, ldr_bot, tau, work )
00752     implicit none
00753 
00754     integer, value, intent(in) :: n, ldr_top, ldr_bot
00755     complex (4), intent(inout) :: R_top(ldr_top,n), R_bot(ldr_bot,n)
00756     complex (4), intent(out)   :: tau(n)
00757     complex (4), intent(out)   :: work(n)
00758 
00759     complex (4)                :: ZERO, ONE
00760     integer                    :: j, k
00761 
00762     ZERO = complex( 0.0e0, 0.0e0 )
00763     ONE = complex( 1.0e0, 0.0e0 )
00764     do k = 1, n
00765        work(k) = ZERO
00766     end do
00767 
00768     do k = 1, n-1
00769        call CLARFP_wrapper( k + 1, R_top(k,k), R_bot(1,k), 1, tau(k) )
00770        call CGEMV( 'Conjugate transpose', k, n-k, ONE, &
00771             R_bot(1, k+1), ldr_bot, &
00772             R_bot(1, k), 1, ZERO, work, 1 )
00773        do j = k+1, n
00774           work(j-k) = work(j-k) + conjg( R_top(k, j) )
00775        end do
00776        do j = k+1, n
00777           R_top(k, j) = R_top(k, j) - tau(k) * work(j-k)
00778        end do
00779        call CGERC( k, n-k, -tau(k), R_bot(1,k), 1, work, 1, &
00780             R_bot(1, k+1), ldr_bot )
00781     end do
00782 
00783     call CLARFP_wrapper( n + 1, R_top(n,n), R_bot(1,n), 1, tau(n) )
00784   end subroutine c_factor_pair
00785 
00786 
00787   ! Apply Q factor (or Q^T) of the 2*ncols_Q by ncols_Q matrix 
00788   ! [R_top; R_bot] (stored in R_bot and tau) to the 2*ncols_Q by ncols_C 
00789   ! matrix [C_top; C_bot].  The two blocks C_top and C_bot may have 
00790   ! different leading dimensions (ldc_top resp. ldc_bot).
00791   !
00792   subroutine d_apply_pair( trans, ncols_C, ncols_Q, R_bot, ldr_bot, &
00793        tau, C_top, ldc_top, C_bot, ldc_bot, work )
00794     implicit none
00795     
00796     character, intent(in)        :: trans
00797     integer, intent(in), value   :: ncols_Q, ncols_C, ldr_bot, ldc_top, ldc_bot
00798     real(8), intent(in)          :: R_bot(ldr_bot,ncols_Q)
00799     real(8), intent(inout)       :: C_top(ldc_top,ncols_C), C_bot(ldc_bot,ncols_C)
00800     real(8), intent(in)          :: tau(ncols_Q)
00801     real(8), intent(out), target :: work(ncols_C)
00802 
00803     real(8), pointer :: y(:)
00804     integer          :: j, k
00805 
00806     y => work(1:ncols_C)
00807     y = 0
00808     if (trans == 'N' .or. trans == 'n') then
00809        do k = ncols_Q, 1, -1
00810           ! The k bottom elements of the length k+1 "sparse" Householder
00811           ! reflector are stored in R_bot(1:k,k).  The topmost element is
00812           ! implicit and is 1.0.  If we want to exploit its sparsity, we
00813           ! can't use DLARF to apply it to [C_top; C_bot]; we have to
00814           ! roll our own routine, which you see here.
00815           
00816           ! $y^T := R_bot(1:k,k)^T C_bot(1:k, 1:ncols_C)$, so $y := C_bot(1:k, 1:ncols_C)^T R_bot(1:k, k)$.
00817           call DGEMV( 'T', k, ncols_C, 1.0d0, C_bot, ldc_bot, R_bot(1, k), 1, 0.0d0, y, 1 )
00818           
00819           ! $y^T := y^T + C_top(k, 1:ncols_C)$
00820           do j = 1, ncols_C
00821              y(j) = y(j) + C_top(k, j)
00822           end do
00823           
00824           ! Update C_top(k, 1:ncols_C)
00825           do j = 1, ncols_C
00826              C_top(k, j) = C_top(k, j) - tau(k) * y(j)
00827           end do
00828           
00829           ! Update C_bot(1:k, 1:ncols_C)
00830           call DGER( k, ncols_C, -tau(k), R_bot(1,k), 1, y, 1, C_bot, ldc_bot )
00831        end do
00832     else
00833        do k = 1, ncols_Q
00834           ! The k bottom elements of the length k+1 "sparse" Householder
00835           ! reflector are stored in R_bot(1:k,k).  The topmost element is
00836           ! implicit and is 1.0.  If we want to exploit its sparsity, we
00837           ! can't use DLARF to apply it to [C_top; C_bot]; we have to
00838           ! roll our own routine, which you see here.
00839           
00840           ! $y^T := R_bot(1:k,k)^T C_bot(1:k, 1:ncols_C)$, so $y := C_bot(1:k, 1:ncols_C)^T R_bot(1:k, k)$.
00841           call DGEMV( 'T', k, ncols_C, 1.0d0, C_bot, ldc_bot, R_bot(1, k), 1, 0.0d0, y, 1 )
00842           
00843           ! $y^T := y^T + C_top(k, 1:ncols_C)$
00844           do j = 1, ncols_C
00845              y(j) = y(j) + C_top(k, j)
00846           end do
00847           
00848           ! Update C_top(k, 1:ncols_C)
00849           do j = 1, ncols_C
00850              C_top(k, j) = C_top(k, j) - tau(k) * y(j)
00851           end do
00852           
00853           ! Update C_bot(1:k, 1:ncols_C)
00854           call DGER( k, ncols_C, -tau(k), R_bot(1,k), 1, y, 1, C_bot(1, 1), ldc_bot )
00855        end do
00856     end if
00857   end subroutine d_apply_pair
00858 
00859 
00860   subroutine z_apply_pair( trans, ncols_C, ncols_Q, R_bot, ldr_bot, &
00861        tau, C_top, ldc_top, C_bot, ldc_bot, work )
00862     implicit none
00863     
00864     character, intent(in)      :: trans
00865     integer, intent(in), value :: ncols_Q, ncols_C, ldr_bot, ldc_top, ldc_bot
00866     complex(8), intent(in)     :: R_bot(ldr_bot,ncols_Q)
00867     complex(8), intent(inout)  :: C_top(ldc_top,ncols_C), C_bot(ldc_bot,ncols_C)
00868     complex(8), intent(in)     :: tau(ncols_Q)
00869     complex(8), intent(out)    :: work(ncols_C)
00870 
00871     complex(8)                 :: ZERO, ONE, tau_k
00872     integer                    :: j, k, k_first, k_second, k_step
00873     logical                    :: no_trans
00874 
00875     ZERO = complex( 0.0d0, 0.0d0 )
00876     ONE = complex( 1.0d0, 0.0d0 )
00877     do k = 1, ncols_C
00878        work(k) = ZERO
00879     end do
00880 
00881     no_trans = (trans == 'N' .or. trans == 'n')
00882     if (no_trans) then
00883        k_first = ncols_Q
00884        k_second = 1
00885        k_step = -1
00886     else 
00887        k_first = 1
00888        k_second = ncols_Q
00889        k_step = 1
00890     end if
00891 
00892     do k = k_first, k_second, k_step
00893        if (no_trans) then
00894           tau_k = tau(k)
00895        else
00896           tau_k = conjg( tau(k) )
00897        end if
00898 
00899        call ZGEMV( 'Conjugate transpose', k, ncols_C, ONE, C_bot, ldc_bot, &
00900             R_bot(1, k), 1, ZERO, work, 1 )
00901        do j = 1, ncols_C
00902           work(j) = work(j) + conjg( C_top(k, j) )
00903        end do
00904        do j = 1, ncols_C
00905           C_top(k, j) = C_top(k, j) - tau_k * work(j)
00906        end do
00907        call ZGERC( k, ncols_C, -tau_k, R_bot(1,k), 1, work, 1, C_bot, ldc_bot )
00908     end do
00909   end subroutine z_apply_pair
00910 
00911 
00912   subroutine s_apply_pair( trans, ncols_C, ncols_Q, R_bot, ldr_bot, &
00913        tau, C_top, ldc_top, C_bot, ldc_bot, work )
00914     implicit none
00915     
00916     character, intent(in)        :: trans
00917     integer, intent(in), value   :: ncols_Q, ncols_C, ldr_bot, ldc_top, ldc_bot
00918     real(4), intent(in)          :: R_bot(ldr_bot,ncols_Q)
00919     real(4), intent(inout)       :: C_top(ldc_top,ncols_C), C_bot(ldc_bot,ncols_C)
00920     real(4), intent(in)          :: tau(ncols_Q)
00921     real(4), intent(out), target :: work(ncols_C)
00922 
00923     real(4), pointer :: y(:)
00924     integer          :: j, k
00925 
00926     y => work(1:ncols_C)
00927     y = 0
00928     if (trans == 'N' .or. trans == 'n') then
00929        do k = ncols_Q, 1, -1
00930           ! The k bottom elements of the length k+1 "sparse" Householder
00931           ! reflector are stored in R_bot(1:k,k).  The topmost element is
00932           ! implicit and is 1.0.  If we want to exploit its sparsity, we
00933           ! can't use DLARF to apply it to [C_top; C_bot]; we have to
00934           ! roll our own routine, which you see here.
00935           
00936           ! $y^T := R_bot(1:k,k)^T C_bot(1:k, 1:ncols_C)$, so $y := C_bot(1:k, 1:ncols_C)^T R_bot(1:k, k)$.
00937           call SGEMV( 'T', k, ncols_C, 1.0, C_bot, ldc_bot, R_bot(1, k), 1, 0.0, y, 1 )
00938           
00939           ! $y^T := y^T + C_top(k, 1:ncols_C)$
00940           do j = 1, ncols_C
00941              y(j) = y(j) + C_top(k, j)
00942           end do
00943           
00944           ! Update C_top(k, 1:ncols_C)
00945           do j = 1, ncols_C
00946              C_top(k, j) = C_top(k, j) - tau(k) * y(j)
00947           end do
00948           
00949           ! Update C_bot(1:k, 1:ncols_C)
00950           call SGER( k, ncols_C, -tau(k), R_bot(1,k), 1, y, 1, C_bot, ldc_bot )
00951        end do
00952     else
00953        do k = 1, ncols_Q
00954           ! The k bottom elements of the length k+1 "sparse" Householder
00955           ! reflector are stored in R_bot(1:k,k).  The topmost element is
00956           ! implicit and is 1.0.  If we want to exploit its sparsity, we
00957           ! can't use DLARF to apply it to [C_top; C_bot]; we have to
00958           ! roll our own routine, which you see here.
00959           
00960           ! $y^T := R_bot(1:k,k)^T C_bot(1:k, 1:ncols_C)$, so $y := C_bot(1:k, 1:ncols_C)^T R_bot(1:k, k)$.
00961           call SGEMV( 'T', k, ncols_C, 1.0, C_bot, ldc_bot, R_bot(1, k), 1, 0.0, y, 1 )
00962           
00963           ! $y^T := y^T + C_top(k, 1:ncols_C)$
00964           do j = 1, ncols_C
00965              y(j) = y(j) + C_top(k, j)
00966           end do
00967           
00968           ! Update C_top(k, 1:ncols_C)
00969           do j = 1, ncols_C
00970              C_top(k, j) = C_top(k, j) - tau(k) * y(j)
00971           end do
00972           
00973           ! Update C_bot(1:k, 1:ncols_C)
00974           call SGER( k, ncols_C, -tau(k), R_bot(1,k), 1, y, 1, C_bot(1, 1), ldc_bot )
00975        end do
00976     end if
00977   end subroutine s_apply_pair
00978 
00979 
00980   subroutine c_apply_pair( trans, ncols_C, ncols_Q, R_bot, ldr_bot, &
00981        tau, C_top, ldc_top, C_bot, ldc_bot, work )
00982     implicit none
00983     
00984     character, intent(in)      :: trans
00985     integer, intent(in), value :: ncols_Q, ncols_C, ldr_bot, ldc_top, ldc_bot
00986     complex(4), intent(in)     :: R_bot(ldr_bot,ncols_Q)
00987     complex(4), intent(inout)  :: C_top(ldc_top,ncols_C), C_bot(ldc_bot,ncols_C)
00988     complex(4), intent(in)     :: tau(ncols_Q)
00989     complex(4), intent(out)    :: work(ncols_C)
00990 
00991     complex(4)                 :: ZERO, ONE, tau_k
00992     integer                    :: j, k, k_first, k_second, k_step
00993     logical                    :: no_trans
00994 
00995     ZERO = complex( 0.0e0, 0.0e0 )
00996     ONE = complex( 1.0e0, 0.0e0 )
00997     do k = 1, ncols_C
00998        work(k) = ZERO
00999     end do
01000 
01001     no_trans = (trans == 'N' .or. trans == 'n')
01002     if (no_trans) then
01003        k_first = ncols_Q
01004        k_second = 1
01005        k_step = -1
01006     else 
01007        k_first = 1
01008        k_second = ncols_Q
01009        k_step = 1
01010     end if
01011 
01012     do k = k_first, k_second, k_step
01013        if (no_trans) then
01014           tau_k = tau(k)
01015        else
01016           tau_k = conjg( tau(k) )
01017        end if
01018 
01019        call CGEMV( 'Conjugate transpose', k, ncols_C, ONE, C_bot, ldc_bot, &
01020             R_bot(1, k), 1, ZERO, work, 1 )
01021        do j = 1, ncols_C
01022           work(j) = work(j) + conjg( C_top(k, j) )
01023        end do
01024        do j = 1, ncols_C
01025           C_top(k, j) = C_top(k, j) - tau_k * work(j)
01026        end do
01027        call CGERC( k, ncols_C, -tau_k, R_bot(1,k), 1, work, 1, C_bot, ldc_bot )
01028     end do
01029   end subroutine c_apply_pair
01030 
01031 
01032 
01033 end module TsqrCombine
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends