ScaLAPACK  2.0.2
ScaLAPACK: Scalable Linear Algebra PACKage
pdtzrzf.f
Go to the documentation of this file.
00001       SUBROUTINE PDTZRZF( M, N, A, IA, JA, DESCA, TAU, WORK, LWORK,
00002      $                    INFO )
00003 *
00004 *  -- ScaLAPACK routine (version 1.7) --
00005 *     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
00006 *     and University of California, Berkeley.
00007 *     May 25, 2001
00008 *
00009 *     .. Scalar Arguments ..
00010       INTEGER            IA, INFO, JA, LWORK, M, N
00011 *     ..
00012 *     .. Array Arguments ..
00013       INTEGER            DESCA( * )
00014       DOUBLE PRECISION   A( * ), TAU( * ), WORK( * )
00015 *     ..
00016 *
00017 *  Purpose
00018 *  =======
00019 *
00020 *  PDTZRZF reduces the M-by-N ( M<=N ) real upper trapezoidal matrix
00021 *  sub( A ) = A(IA:IA+M-1,JA:JA+N-1) to upper triangular form by means
00022 *  of orthogonal transformations.
00023 *
00024 *  The upper trapezoidal matrix sub( A ) is factored as
00025 *
00026 *     sub( A ) = ( R  0 ) * Z,
00027 *
00028 *  where Z is an N-by-N orthogonal matrix and R is an M-by-M upper
00029 *  triangular matrix.
00030 *
00031 *  Notes
00032 *  =====
00033 *
00034 *  Each global data object is described by an associated description
00035 *  vector.  This vector stores the information required to establish
00036 *  the mapping between an object element and its corresponding process
00037 *  and memory location.
00038 *
00039 *  Let A be a generic term for any 2D block cyclicly distributed array.
00040 *  Such a global array has an associated description vector DESCA.
00041 *  In the following comments, the character _ should be read as
00042 *  "of the global array".
00043 *
00044 *  NOTATION        STORED IN      EXPLANATION
00045 *  --------------- -------------- --------------------------------------
00046 *  DTYPE_A(global) DESCA( DTYPE_ )The descriptor type.  In this case,
00047 *                                 DTYPE_A = 1.
00048 *  CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating
00049 *                                 the BLACS process grid A is distribu-
00050 *                                 ted over. The context itself is glo-
00051 *                                 bal, but the handle (the integer
00052 *                                 value) may vary.
00053 *  M_A    (global) DESCA( M_ )    The number of rows in the global
00054 *                                 array A.
00055 *  N_A    (global) DESCA( N_ )    The number of columns in the global
00056 *                                 array A.
00057 *  MB_A   (global) DESCA( MB_ )   The blocking factor used to distribute
00058 *                                 the rows of the array.
00059 *  NB_A   (global) DESCA( NB_ )   The blocking factor used to distribute
00060 *                                 the columns of the array.
00061 *  RSRC_A (global) DESCA( RSRC_ ) The process row over which the first
00062 *                                 row of the array A is distributed.
00063 *  CSRC_A (global) DESCA( CSRC_ ) The process column over which the
00064 *                                 first column of the array A is
00065 *                                 distributed.
00066 *  LLD_A  (local)  DESCA( LLD_ )  The leading dimension of the local
00067 *                                 array.  LLD_A >= MAX(1,LOCr(M_A)).
00068 *
00069 *  Let K be the number of rows or columns of a distributed matrix,
00070 *  and assume that its process grid has dimension p x q.
00071 *  LOCr( K ) denotes the number of elements of K that a process
00072 *  would receive if K were distributed over the p processes of its
00073 *  process column.
00074 *  Similarly, LOCc( K ) denotes the number of elements of K that a
00075 *  process would receive if K were distributed over the q processes of
00076 *  its process row.
00077 *  The values of LOCr() and LOCc() may be determined via a call to the
00078 *  ScaLAPACK tool function, NUMROC:
00079 *          LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ),
00080 *          LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ).
00081 *  An upper bound for these quantities may be computed by:
00082 *          LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A
00083 *          LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A
00084 *
00085 *  Arguments
00086 *  =========
00087 *
00088 *  M       (global input) INTEGER
00089 *          The number of rows to be operated on, i.e. the number of rows
00090 *          of the distributed submatrix sub( A ). M >= 0.
00091 *
00092 *  N       (global input) INTEGER
00093 *          The number of columns to be operated on, i.e. the number of
00094 *          columns of the distributed submatrix sub( A ). N >= 0.
00095 *
00096 *  A       (local input/local output) DOUBLE PRECISION pointer into the
00097 *          local memory to an array of dimension (LLD_A, LOCc(JA+N-1)).
00098 *          On entry, the local pieces of the M-by-N distributed matrix
00099 *          sub( A ) which is to be factored. On exit, the leading M-by-M
00100 *          upper triangular part of sub( A ) contains the upper trian-
00101 *          gular matrix R, and elements M+1 to N of the first M rows of
00102 *          sub( A ), with the array TAU, represent the orthogonal matrix
00103 *          Z as a product of M elementary reflectors.
00104 *
00105 *  IA      (global input) INTEGER
00106 *          The row index in the global array A indicating the first
00107 *          row of sub( A ).
00108 *
00109 *  JA      (global input) INTEGER
00110 *          The column index in the global array A indicating the
00111 *          first column of sub( A ).
00112 *
00113 *  DESCA   (global and local input) INTEGER array of dimension DLEN_.
00114 *          The array descriptor for the distributed matrix A.
00115 *
00116 *  TAU     (local output) DOUBLE PRECISION array, dimension LOCr(IA+M-1)
00117 *          This array contains the scalar factors of the elementary
00118 *          reflectors. TAU is tied to the distributed matrix A.
00119 *
00120 *  WORK    (local workspace/local output) DOUBLE PRECISION array,
00121 *                                                    dimension (LWORK)
00122 *          On exit, WORK(1) returns the minimal and optimal LWORK.
00123 *
00124 *  LWORK   (local or global input) INTEGER
00125 *          The dimension of the array WORK.
00126 *          LWORK is local input and must be at least
00127 *          LWORK >= MB_A * ( Mp0 + Nq0 + MB_A ), where
00128 *
00129 *          IROFF = MOD( IA-1, MB_A ), ICOFF = MOD( JA-1, NB_A ),
00130 *          IAROW = INDXG2P( IA, MB_A, MYROW, RSRC_A, NPROW ),
00131 *          IACOL = INDXG2P( JA, NB_A, MYCOL, CSRC_A, NPCOL ),
00132 *          Mp0   = NUMROC( M+IROFF, MB_A, MYROW, IAROW, NPROW ),
00133 *          Nq0   = NUMROC( N+ICOFF, NB_A, MYCOL, IACOL, NPCOL ),
00134 *
00135 *          and NUMROC, INDXG2P are ScaLAPACK tool functions;
00136 *          MYROW, MYCOL, NPROW and NPCOL can be determined by calling
00137 *          the subroutine BLACS_GRIDINFO.
00138 *
00139 *          If LWORK = -1, then LWORK is global input and a workspace
00140 *          query is assumed; the routine only calculates the minimum
00141 *          and optimal size for all work arrays. Each of these
00142 *          values is returned in the first entry of the corresponding
00143 *          work array, and no error message is issued by PXERBLA.
00144 *
00145 *  INFO    (global output) INTEGER
00146 *          = 0:  successful exit
00147 *          < 0:  If the i-th argument is an array and the j-entry had
00148 *                an illegal value, then INFO = -(i*100+j), if the i-th
00149 *                argument is a scalar and had an illegal value, then
00150 *                INFO = -i.
00151 *
00152 *  Further Details
00153 *  ===============
00154 *
00155 *  The  factorization is obtained by Householder's method.  The kth
00156 *  transformation matrix, Z( k ), which is used to introduce zeros into
00157 *  the (m - k + 1)th row of sub( A ), is given in the form
00158 *
00159 *     Z( k ) = ( I     0   ),
00160 *              ( 0  T( k ) )
00161 *
00162 *  where
00163 *
00164 *     T( k ) = I - tau*u( k )*u( k )',   u( k ) = (   1    ),
00165 *                                                 (   0    )
00166 *                                                 ( z( k ) )
00167 *
00168 *  tau is a scalar and z( k ) is an ( n - m ) element vector.
00169 *  tau and z( k ) are chosen to annihilate the elements of the kth row
00170 *  of sub( A ).
00171 *
00172 *  The scalar tau is returned in the kth element of TAU and the vector
00173 *  u( k ) in the kth row of sub( A ), such that the elements of z( k )
00174 *  are in  a( k, m + 1 ), ..., a( k, n ). The elements of R are returned
00175 *  in the upper triangular part of sub( A ).
00176 *
00177 *  Z is given by
00178 *
00179 *     Z =  Z( 1 ) * Z( 2 ) * ... * Z( m ).
00180 *
00181 *  =====================================================================
00182 *
00183 *     .. Parameters ..
00184       INTEGER            BLOCK_CYCLIC_2D, CSRC_, CTXT_, DLEN_, DTYPE_,
00185      $                   LLD_, MB_, M_, NB_, N_, RSRC_
00186       PARAMETER          ( BLOCK_CYCLIC_2D = 1, DLEN_ = 9, DTYPE_ = 1,
00187      $                     CTXT_ = 2, M_ = 3, N_ = 4, MB_ = 5, NB_ = 6,
00188      $                     RSRC_ = 7, CSRC_ = 8, LLD_ = 9 )
00189       DOUBLE PRECISION   ZERO
00190       PARAMETER          ( ZERO = 0.0D+0 )
00191 *     ..
00192 *     .. Local Scalars ..
00193       LOGICAL            LQUERY
00194       CHARACTER          COLBTOP, ROWBTOP
00195       INTEGER            I, IACOL, IAROW, IB, ICTXT, IIA, IL, IN, IPW,
00196      $                   IROFFA, J, JM1, L, LWMIN, MP0, MYCOL, MYROW,
00197      $                   NPCOL, NPROW, NQ0
00198 *     ..
00199 *     .. Local Arrays ..
00200       INTEGER            IDUM1( 1 ), IDUM2( 1 )
00201 *     ..
00202 *     .. External Subroutines ..
00203       EXTERNAL           BLACS_GRIDINFO, CHK1MAT, INFOG1L, PCHK1MAT,
00204      $                   PDLATRZ, PDLARZB, PDLARZT, PB_TOPGET,
00205      $                   PB_TOPSET, PXERBLA
00206 *     ..
00207 *     .. External Functions ..
00208       INTEGER            ICEIL, INDXG2P, NUMROC
00209       EXTERNAL           ICEIL, INDXG2P, NUMROC
00210 *     ..
00211 *     .. Intrinsic Functions ..
00212       INTRINSIC          DBLE, MAX, MIN, MOD
00213 *     ..
00214 *     .. Executable Statements ..
00215 *
00216 *     Get grid parameters
00217 *
00218       ICTXT = DESCA( CTXT_ )
00219       CALL BLACS_GRIDINFO( ICTXT, NPROW, NPCOL, MYROW, MYCOL )
00220 *
00221 *     Test the input parameters
00222 *
00223       INFO = 0
00224       IF( NPROW.EQ.-1 ) THEN
00225          INFO = -(600+CTXT_)
00226       ELSE
00227          CALL CHK1MAT( M, 1, N, 2, IA, JA, DESCA, 6, INFO )
00228          IF( INFO.EQ.0 ) THEN
00229             IROFFA = MOD( IA-1, DESCA( MB_ ) )
00230             IAROW = INDXG2P( IA, DESCA( MB_ ), MYROW, DESCA( RSRC_ ),
00231      $                       NPROW )
00232             IACOL = INDXG2P( JA, DESCA( NB_ ), MYCOL, DESCA( CSRC_ ),
00233      $                       NPCOL )
00234             MP0 = NUMROC( M+IROFFA, DESCA( MB_ ), MYROW, IAROW, NPROW )
00235             NQ0 = NUMROC( N+MOD( JA-1, DESCA( NB_ ) ), DESCA( NB_ ),
00236      $                    MYCOL, IACOL, NPCOL )
00237             LWMIN = DESCA( MB_ ) * ( MP0 + NQ0 + DESCA( MB_ ) )
00238 *
00239             WORK( 1 ) = DBLE( LWMIN )
00240             LQUERY = ( LWORK.EQ.-1 )
00241             IF( N.LT.M ) THEN
00242                INFO = -2
00243             ELSE IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
00244                INFO = -9
00245             END IF
00246          END IF
00247          IF( LQUERY ) THEN
00248             IDUM1( 1 ) = -1
00249          ELSE
00250             IDUM1( 1 ) = 1
00251          END IF
00252          IDUM2( 1 ) = 9
00253          CALL PCHK1MAT( M, 1, N, 2, IA, JA, DESCA, 6, 1, IDUM1, IDUM2,
00254      $                  INFO )
00255       END IF
00256 *
00257       IF( INFO.NE.0 ) THEN
00258          CALL PXERBLA( ICTXT, 'PDTZRZF', -INFO )
00259          RETURN
00260       ELSE IF( LQUERY ) THEN
00261          RETURN
00262       END IF
00263 *
00264 *     Quick return if possible
00265 *
00266       IF( M.EQ.0 .OR. N.EQ.0 )
00267      $   RETURN
00268 *
00269       IF( M.EQ.N ) THEN
00270 *
00271          CALL INFOG1L( IA, DESCA( MB_ ), NPROW, MYROW, DESCA( RSRC_ ),
00272      $                 IIA, IAROW )
00273          IF( MYROW.EQ.IAROW )
00274      $      MP0 = MP0 - IROFFA
00275          DO 10 I = IIA, IIA+MP0-1
00276             TAU( I ) = ZERO
00277    10    CONTINUE
00278 *
00279       ELSE
00280 *
00281          L = N-M
00282          JM1 = JA + MIN( M+1, N ) - 1
00283          IPW = DESCA( MB_ ) * DESCA( MB_ ) + 1
00284          IN = MIN( ICEIL( IA, DESCA( MB_ ) ) * DESCA( MB_ ), IA+M-1 )
00285          IL = MAX( ( (IA+M-2) / DESCA( MB_ ) ) * DESCA( MB_ ) + 1, IA )
00286          CALL PB_TOPGET( ICTXT, 'Broadcast', 'Rowwise', ROWBTOP )
00287          CALL PB_TOPGET( ICTXT, 'Broadcast', 'Columnwise', COLBTOP )
00288          CALL PB_TOPSET( ICTXT, 'Broadcast', 'Rowwise', ' ' )
00289          CALL PB_TOPSET( ICTXT, 'Broadcast', 'Columnwise', 'D-ring' )
00290 *
00291 *        Use blocked code initially
00292 *
00293          DO 20 I = IL, IN+1, -DESCA( MB_ )
00294             IB = MIN( IA+M-I, DESCA( MB_ ) )
00295             J = JA + I - IA
00296 *
00297 *           Compute the complete orthogonal factorization of the current
00298 *           block A(i:i+ib-1,j:ja+n-1)
00299 *
00300             CALL PDLATRZ( IB, JA+N-J, L, A, I, J, DESCA, TAU, WORK )
00301 *
00302             IF( I.GT.IA ) THEN
00303 *
00304 *              Form the triangular factor of the block reflector
00305 *              H = H(i+ib-1) . . . H(i+1) H(i)
00306 *
00307                CALL PDLARZT( 'Backward', 'Rowwise', L, IB, A, I, JM1,
00308      $                       DESCA, TAU, WORK, WORK( IPW ) )
00309 *
00310 *              Apply H to A(ia:i-1,j:ja+n-1) from the right
00311 *
00312                CALL PDLARZB( 'Right', 'No transpose', 'Backward',
00313      $                       'Rowwise', I-IA, JA+N-J, IB, L, A, I, JM1,
00314      $                       DESCA, WORK, A, IA, J, DESCA, WORK( IPW ) )
00315             END IF
00316 *
00317    20    CONTINUE
00318 *
00319 *        Use unblocked code to factor the last or only block
00320 *
00321          CALL PDLATRZ( IN-IA+1, N, N-M, A, IA, JA, DESCA, TAU, WORK )
00322 *
00323          CALL PB_TOPSET( ICTXT, 'Broadcast', 'Rowwise', ROWBTOP )
00324          CALL PB_TOPSET( ICTXT, 'Broadcast', 'Columnwise', COLBTOP )
00325 *
00326       END IF
00327 *
00328       WORK( 1 ) = DBLE( LWMIN )
00329 *
00330       RETURN
00331 *
00332 *     End of PDTZRZF
00333 *
00334       END