Routine: PCDBTRF()  File: SRC\pcdbtrf.f

 
 
# lines: 1268
  # code: 1268
  # comment: 0
  # blank:0
# Variables:32
# Callers:1
# Callings:0
# Words:300
# Keywords:70
 

 

 
..
     .. Local Scalars ..
     ..
     .. Local Arrays ..
     ..
     .. External Subroutines ..
     ..
     .. External Functions ..
     ..
     .. Intrinsic Functions ..
     ..
     .. Executable Statements ..
     Test the input parameters
     Convert descriptor into standard form for easy access to
        parameters, check that grid is of right shape.
     Get values out of descriptor for use in code.
     Get grid parameters
     Size of separator blocks is maximum of bandwidths
     Argument checking that is specific to Divide & Conquer routine
     Check auxiliary storage size
        put minimum value of laf into AF( 1 )
     Check worksize
     Pack params and positions into arrays for global consistency check
     Want to find errors with MIN( ), so if no error, set it to a big
     number. If there already is an error, multiply by the the
     descriptor multiplier.
     Check consistency across processors
     Prepare output: set info = 0 if no error, and divide by DESCMULT
     if error is not in a descriptor entry.
     Quick return if possible
     Adjust addressing into matrix space to properly get into
        the beginning part of the relevant data
     Form a new BLACS grid (the "standard form" grid) with only procs
        holding part of the matrix, of size 1xNP where NP is adjusted,
        starting at csrc=0, with JA modified to reflect dropped procs.
     First processor to hold part of the matrix:
     Calculate new JA one while dropping off unused processors.
     Save and compute new value of NP
     Call utility routine that forms "standard-form" grid
     Use new context from standard grid as context.
     Get information about new grid.
     Drop out processors that do not have part of the matrix.
     ********************************
     Values reused throughout routine
     User-input value of partition size
     Number of columns in each processor
     Offset in columns to beginning of main partition in each proc
     Offset in elements
     Size of main (or odd) partition in each processor
     Offset to workspace for Upper triangular factor
       Zero out space for fillin
       Zero out space for work
     Begin main code
*******************************************************************
       PHASE 1: Local computation phase.
*******************************************************************
       Sizes of the extra triangles communicated bewtween processors
         Transfer last triangle D_i of local matrix to next processor
         which needs it to calculate fillin due to factorization of
         its main (odd) block A_i.
         Overlap the send with the factorization of A_i.
       Factor main partition A_i = L_i {U_i} in each processor
         Apply factorization to lower connection block BL_i
         conjugate transpose the connection block in preparation.
         Apply factorization to upper connection block BU_i
         Move the connection block in preparation.
         Perform the triangular system solve {L_i}{{BU'}_i} = {B_i}
         Perform the triangular solve {U_i}^C{BL'}_i^C = {BL_i}^C
         conjugate transpose resulting block to its location
           in main storage.
         Move the resulting block back to its location in main storage.
         Compute contribution to diagonal block(s) of reduced system.
          {C'}_i = {C_i}-{{BL'}_i}{{BU'}_i}
         The following method uses more flops than necessary but
           does not necessitate the writing of a new BLAS routine.
       End of "if ( MYCOL .lt. NP-1 )..." loop
       If the processor could not locally factor, it jumps here.
         Discard temporary matrix stored beginning in
           AF( (odd_size+2*bwl, bwu)*bwl, bwu+1 ) and use for
           off_diagonal block of reduced system.
         Receive previously transmitted matrix section, which forms
         the right-hand-side for the triangular solve that calculates
         the "spike" fillin.
         Calculate the "spike" fillin, ${L_i} {{GU}_i} = {DL_i}$ .
         Calculate the "spike" fillin, ${U_i}^C {{GL}_i}^C = {DU_i}^C$
         Copy D block into AF storage for solve.
         Calculate the update block for previous proc, E_i = GL_i{GU_i}
         Zero out space in case result is smaller than storage block
         Initiate send of E_i to previous processor to overlap
           with next computation.
           Calculate off-diagonal block(s) of reduced system.
           Note: for ease of use in solution of reduced system, store
           L's off-diagonal block in conjugate transpose form.
           Copy matrix HU_i (the last bwl rows of GU_i) to AFL storage
             as per requirements of BLAS routine CTRMM.
             Since we have GU_i stored,
             conjugate transpose HU_i to HU_i^C.
           Copy matrix HL_i (the last bwu rows of GL_i^C) to AFU store
             as per requirements of BLAS routine CTRMM.
             Since we have GL_i^C stored,
             conjugate transpose HL_i^C to HL_i.
       End of "if ( MYCOL .ne. 0 )..."
       End of "if (info.eq.0) then"
       Check to make sure no processors have found errors
       No errors found, continue
*******************************************************************
       PHASE 2: Formation and factorization of Reduced System.
*******************************************************************
       Gather up local sections of reduced system
     The last processor does not participate in the factorization of
       the reduced system, having sent its E_i already.
       Initiate send of off-diag block(s) to overlap with next part.
       Off-diagonal block needed on neighboring processor to start
       algorithm.
       Copy last diagonal block into AF storage for subsequent
         operations.
       Receive cont. to diagonal block that is stored on this proc.
          Add contribution to diagonal block
       *************************************
       Modification Loop
       The distance for sending and receiving for each level starts
         at 1 for the first level.
       Do until this proc is needed to modify other procs' equations
         Receive and add contribution to diagonal block from the left
         Receive and add contribution to diagonal block from the right
       [End of GOTO Loop]
       *********************************
       Calculate and use this proc's blocks to modify other procs'...
       Factor diagonal block
       ****************************************************************
       Receive offdiagonal block from processor to right.
         If this is the first group of processors, the receive comes
         from a different processor than otherwise.
           Move block into place that it will be expected to be for
             calcs.
 
Display dynamic version Find AutoScroll Reload FontSize: - + Hide Comments Hide Blanks Frame FullScreen MailPrint

  
001        SUBROUTINE PCDBTRF( N , BWL , BWU , A , JA , DESCA , AF , LAF , WORK ,
002       $LWORK , INFO )
003  
004  *     -- ScaLAPACK routine(version 1.7) --
005  *     University of Tennessee , Knoxville , Oak Ridge National Laboratory ,
006  *     and University of California , Berkeley.
007  *     May 25 , 2001
008  
009  *     .. Scalar Arguments ..
010        INTEGER BWL , BWU , INFO , JA , LAF , LWORK , N
011  *     ..
012  *     .. Array Arguments ..
013        INTEGER DESCA( * )
014        COMPLEX A( * ) , AF( * ) , WORK( * )
015  *     ..
016  
017  *     Purpose
018  *     === ====
019  
020  *     PCDBTRF computes a LU factorization
021  *     of an N - by - N complex banded
022  *     diagonally dominant - like distributed matrix
023  *     with bandwidth BWL , BWU : A(1 : N , JA : JA + N - 1).
024  *     Reordering is used to increase parallelism in the factorization.
025  *     This reordering results in factors that are DIFFERENT from those
026  *     produced by equivalent sequential codes. These factors cannot
027  *     be used directly by users ; however , they can be used in
028  *     subsequent calls to PCDBTRS to solve linear systems.
029  
030  *     The factorization has the form
031  
032  *     P A(1 : N , JA : JA + N - 1) P^T = L U
033  
034  *     where U is a banded upper triangular matrix and L is banded
035  *     lower triangular , and P is a permutation matrix.
036  
037  *     === ==================================================================
038  
039  *     Arguments
040  *     === ======
041  
042  *     N(global input) INTEGER
043  *     The number of rows and columns to be operated on , i.e. the
044  *     order of the distributed submatrix A(1 : N , JA : JA + N - 1). N >= 0.
045  
046  *     BWL(global input) INTEGER
047  *     Number of subdiagonals. 0 <= BWL <= N - 1
048  
049  *     BWU(global input) INTEGER
050  *     Number of superdiagonals. 0 <= BWU <= N - 1
051  
052  *     A(local input / local output) COMPLEX pointer into
053  *     local memory to an array with first dimension
054  *     LLD_A >=(bwl + bwu + 1)(stored in DESCA).
055  *     On entry , this array contains the local pieces of the
056  *     N - by - N unsymmetric banded distributed matrix
057  *     A(1 : N , JA : JA + N - 1) to be factored.
058  *     This local portion is stored in the packed banded format
059  *     used in LAPACK. Please see the Notes below and the
060  *     ScaLAPACK manual for more detail on the format of
061  *     distributed matrices.
062  *     On exit , this array contains information containing details
063  *     of the factorization.
064  *     Note that permutations are performed on the matrix , so that
065  *     the factors returned are different from those returned
066  *     by LAPACK.
067  
068  *     JA(global input) INTEGER
069  *     The index in the global array A that points to the start of
070  *     the matrix to be operated on(which may be either all of A
071  *     or a submatrix of A).
072  
073  *     DESCA(global and local input) INTEGER array of dimension DLEN.
074  *     if 1D type(DTYPE_A = 501) , DLEN >= 7 ;
075  *     if 2D type(DTYPE_A = 1) , DLEN >= 9 .
076  *     The array descriptor for the distributed matrix A.
077  *     Contains information of mapping of A to memory. Please
078  *     see NOTES below for full description and options.
079  
080  *     AF(local output) COMPLEX array , dimension LAF.
081  *     Auxiliary Fillin Space.
082  *     Fillin is created during the factorization routine
083  *     PCDBTRF and this is stored in AF. If a linear system
084  *     is to be solved using PCDBTRS after the factorization
085  *     routine , AF *must not be altered* after the factorization.
086  
087  *     LAF(local input) INTEGER
088  *     Size of user - input Auxiliary Fillin space AF. Must be >=
089  *     NB*(bwl + bwu) + 6*max(bwl , bwu)*max(bwl , bwu)
090  *     If LAF is not large enough , an error code will be returned
091  *     and the minimum acceptable size will be returned in AF( 1 )
092  
093  *     WORK(local workspace / local output)
094  *     COMPLEX temporary workspace. This space may
095  *     be overwritten in between calls to routines. WORK must be
096  *     the size given in LWORK.
097  *     On exit , WORK( 1 ) contains the minimal LWORK.
098  
099  *     LWORK(local input or global input) INTEGER
100  *     Size of user - input workspace WORK.
101  *     If LWORK is too small , the minimal acceptable size will be
102  *     returned in WORK(1) and an error code is returned. LWORK >=
103  *     max(bwl , bwu)*max(bwl , bwu)
104  
105  *     INFO(global output) INTEGER
106  *     = 0 : successful exit
107  *     < 0 : If the i - th argument is an array and the j - entry had
108  *     an illegal value , then INFO = - (i*100 + j) , if the i - th
109  *     argument is a scalar and had an illegal value , then
110  *     INFO = - i.
111  *     > 0 : If INFO = K <= NPROCS , the submatrix stored on processor
112  *     INFO and factored locally was not
113  *     diagonally dominant - like , and
114  *     the factorization was not completed.
115  *     If INFO = K > NPROCS , the submatrix stored on processor
116  *     INFO - NPROCS representing interactions with other
117  *     processors was not
118  *     stably factorable wo / interchanges ,
119  *     and the factorization was not completed.
120  
121  *     === ==================================================================
122  
123  *     Restrictions
124  *     === =========
125  
126  *     The following are restrictions on the input parameters. Some of these
127  *     are temporary and will be removed in future releases , while others
128  *     may reflect fundamental technical limitations.
129  
130  *     Non - cyclic restriction : VERY IMPORTANT !
131  *     P*NB >= mod(JA - 1 , NB) + N.
132  *     The mapping for matrices must be blocked , reflecting the nature
133  *     of the divide and conquer algorithm as a task - parallel algorithm.
134  *     This formula in words is : no processor may have more than one
135  *     chunk of the matrix.
136  
137  *     Blocksize cannot be too small :
138  *     If the matrix spans more than one processor , the following
139  *     restriction on NB , the size of each block on each processor ,
140  *     must hold :
141  *     NB >= 2*MAX(BWL , BWU)
142  *     The bulk of parallel computation is done on the matrix of size
143  *     O(NB) on each processor. If this is too small , divide and conquer
144  *     is a poor choice of algorithm.
145  
146  *     Submatrix reference :
147  *     JA = IB
148  *     Alignment restriction that prevents unnecessary communication.
149  
150  *     === ==================================================================
151  
152  *     Notes
153  *     === ==
154  
155  *     If the factorization routine and the solve routine are to be called
156  *     separately(to solve various sets of righthand sides using the same
157  *     coefficient matrix) , the auxiliary space AF *must not be altered*
158  *     between calls to the factorization routine and the solve routine.
159  
160  *     The best algorithm for solving banded and tridiagonal linear systems
161  *     depends on a variety of parameters , especially the bandwidth.
162  *     Currently , only algorithms designed for the case N / P >> bw are
163  *     implemented. These go by many names , including Divide and Conquer ,
164  *     Partitioning , domain decomposition - type , etc.
165  
166  *     Algorithm description : Divide and Conquer
167  
168  *     The Divide and Conqer algorithm assumes the matrix is narrowly
169  *     banded compared with the number of equations. In this situation ,
170  *     it is best to distribute the input matrix A one - dimensionally ,
171  *     with columns atomic and rows divided amongst the processes.
172  *     The basic algorithm divides the banded matrix up into
173  *     P pieces with one stored on each processor ,
174  *     and then proceeds in 2 phases for the factorization or 3 for the
175  *     solution of a linear system.
176  *     1) Local Phase :
177  *     The individual pieces are factored independently and in
178  *     parallel. These factors are applied to the matrix creating
179  *     fillin , which is stored in a non - inspectable way in auxiliary
180  *     space AF. Mathematically , this is equivalent to reordering
181  *     the matrix A as P A P^T and then factoring the principal
182  *     leading submatrix of size equal to the sum of the sizes of
183  *     the matrices factored on each processor. The factors of
184  *     these submatrices overwrite the corresponding parts of A
185  *     in memory.
186  *     2) Reduced System Phase :
187  *     A small(max(bwl , bwu)* (P - 1)) system is formed representing
188  *     interaction of the larger blocks , and is stored(as are its
189  *     factors) in the space AF. A parallel Block Cyclic Reduction
190  *     algorithm is used. For a linear system , a parallel front solve
191  *     followed by an analagous backsolve , both using the structure
192  *     of the factored matrix , are performed.
193  *     3) Backsubsitution Phase :
194  *     For a linear system , a local backsubstitution is performed on
195  *     each processor in parallel.
196  
197  *     Descriptors
198  *     === ========
199  
200  *     Descriptors now have *types* and differ from ScaLAPACK 1.0.
201  
202  *     Note : banded codes can use either the old two dimensional
203  *     or new one - dimensional descriptors , though the processor grid in
204  *     both cases *must be one - dimensional*. We describe both types below.
205  
206  *     Each global data object is described by an associated description
207  *     vector. This vector stores the information required to establish
208  *     the mapping between an object element and its corresponding process
209  *     and memory location.
210  
211  *     Let A be a generic term for any 2D block cyclicly distributed array.
212  *     Such a global array has an associated description vector DESCA.
213  *     In the following comments , the character _ should be read as
214  *     "of the global array".
215  
216  *     NOTATION STORED IN EXPLANATION
217  *     --- ------------ -------------- --------------------------------------
218  *     DTYPE_A(global) DESCA( DTYPE_ )The descriptor type. In this case ,
219  *     DTYPE_A = 1.
220  *     CTXT_A(global) DESCA( CTXT_ ) The BLACS context handle , indicating
221  *     the BLACS process grid A is distribu -
222  *     ted over. The context itself is glo -
223  *     bal , but the handle(the integer
224  *     value) may vary.
225  *     M_A(global) DESCA( M_ ) The number of rows in the global
226  *     array A.
227  *     N_A(global) DESCA( N_ ) The number of columns in the global
228  *     array A.
229  *     MB_A(global) DESCA( MB_ ) The blocking factor used to distribute
230  *     the rows of the array.
231  *     NB_A(global) DESCA( NB_ ) The blocking factor used to distribute
232  *     the columns of the array.
233  *     RSRC_A(global) DESCA( RSRC_ ) The process row over which the first
234  *     row of the array A is distributed.
235  *     CSRC_A(global) DESCA( CSRC_ ) The process column over which the
236  *     first column of the array A is
237  *     distributed.
238  *     LLD_A(local) DESCA( LLD_ ) The leading dimension of the local
239  *     array. LLD_A >= MAX(1 , LOCr(M_A)).
240  
241  *     Let K be the number of rows or columns of a distributed matrix ,
242  *     and assume that its process grid has dimension p x q.
243  *     LOCr( K ) denotes the number of elements of K that a process
244  *     would receive if K were distributed over the p processes of its
245  *     process column.
246  *     Similarly , LOCc( K ) denotes the number of elements of K that a
247  *     process would receive if K were distributed over the q processes of
248  *     its process row.
249  *     The values of LOCr() and LOCc() may be determined via a call to the
250  *     ScaLAPACK tool function , NUMROC :
251  *     LOCr( M ) = NUMROC( M , MB_A , MYROW , RSRC_A , NPROW ) ,
252  *     LOCc( N ) = NUMROC( N , NB_A , MYCOL , CSRC_A , NPCOL ).
253  *     An upper bound for these quantities may be computed by :
254  *     LOCr( M ) <= ceil( ceil(M / MB_A) / NPROW )*MB_A
255  *     LOCc( N ) <= ceil( ceil(N / NB_A) / NPCOL )*NB_A
256  
257  *     One - dimensional descriptors :
258  
259  *     One - dimensional descriptors are a new addition to ScaLAPACK since
260  *     version 1.0. They simplify and shorten the descriptor for 1D
261  *     arrays.
262  
263  *     Since ScaLAPACK supports two - dimensional arrays as the fundamental
264  *     object , we allow 1D arrays to be distributed either over the
265  *     first dimension of the array(as if the grid were P - by - 1) or the
266  *     2nd dimension(as if the grid were 1 - by - P). This choice is
267  *     indicated by the descriptor type(501 or 502)
268  *     as described below.
269  
270  *     IMPORTANT NOTE : the actual BLACS grid represented by the
271  *     CTXT entry in the descriptor may be *either* P - by - 1 or 1 - by - P
272  *     irrespective of which one - dimensional descriptor type
273  *     (501 or 502) is input.
274  *     This routine will interpret the grid properly either way.
275  *     ScaLAPACK routines *do not support intercontext operations* so that
276  *     the grid passed to a single ScaLAPACK routine *must be the same*
277  *     for all array descriptors passed to that routine.
278  
279  *     NOTE : In all cases where 1D descriptors are used , 2D descriptors
280  *     may also be used , since a one - dimensional array is a special case
281  *     of a two - dimensional array with one dimension of size unity.
282  *     The two - dimensional array used in this case *must* be of the
283  *     proper orientation :
284  *     If the appropriate one - dimensional descriptor is DTYPEA = 501
285  *     (1 by P type) , then the two dimensional descriptor must
286  *     have a CTXT value that refers to a 1 by P BLACS grid ;
287  *     If the appropriate one - dimensional descriptor is DTYPEA = 502
288  *     (P by 1 type) , then the two dimensional descriptor must
289  *     have a CTXT value that refers to a P by 1 BLACS grid.
290  
291  *     Summary of allowed descriptors , types , and BLACS grids :
292  *     DTYPE 501 502 1 1
293  *     BLACS grid 1xP or Px1 1xP or Px1 1xP Px1
294  *     --- --------------------------------------------------
295  *     A               OK NO OK NO
296  *     B               NO OK NO OK
297  
298  *     Let A be a generic term for any 1D block cyclicly distributed array.
299  *     Such a global array has an associated description vector DESCA.
300  *     In the following comments , the character _ should be read as
301  *     "of the global array".
302  
303  *     NOTATION STORED IN EXPLANATION
304  *     --- ------------ ---------- ------------------------------------------
305  *     DTYPE_A(global) DESCA( 1 ) The descriptor type. For 1D grids ,
306  *     TYPE_A = 501 : 1 - by - P grid.
307  *     TYPE_A = 502 : P - by - 1 grid.
308  *     CTXT_A(global) DESCA( 2 ) The BLACS context handle , indicating
309  *     the BLACS process grid A is distribu -
310  *     ted over. The context itself is glo -
311  *     bal , but the handle(the integer
312  *     value) may vary.
313  *     N_A(global) DESCA( 3 ) The size of the array dimension being
314  *     distributed.
315  *     NB_A(global) DESCA( 4 ) The blocking factor used to distribute
316  *     the distributed dimension of the array.
317  *     SRC_A(global) DESCA( 5 ) The process row or column over which the
318  *     first row or column of the array
319  *     is distributed.
320  *     LLD_A(local) DESCA( 6 ) The leading dimension of the local array
321  *     storing the local blocks of the distri -
322  *     buted array A. Minimum value of LLD_A
323  *     depends on TYPE_A.
324  *     TYPE_A = 501 : LLD_A >=
325  *     size of undistributed dimension , 1.
326  *     TYPE_A = 502 : LLD_A >= NB_A , 1.
327  *     Reserved DESCA( 7 ) Reserved for future use.
328  
329  *     === ==================================================================
330  
331  *     Code Developer : Andrew J. Cleary , University of Tennessee.
332  *     Current address : Lawrence Livermore National Labs.
333  *     This version released : August , 2001.
334  
335  *     === ==================================================================
336  
337  *     ..
338  *     .. Parameters ..
339        REAL ONE , ZERO
340        PARAMETER( ONE = 1.0E + 0 )
341        PARAMETER( ZERO = 0.0E + 0 )
342        COMPLEX CONE , CZERO
343        PARAMETER( CONE =( 1.0E + 0 , 0.0E + 0 ) )
344        PARAMETER( CZERO =( 0.0E + 0 , 0.0E + 0 ) )
345        INTEGER INT_ONE
346        PARAMETER( INT_ONE = 1 )
347        INTEGER DESCMULT , BIGNUM
348        PARAMETER(DESCMULT = 100 , BIGNUM = DESCMULT * DESCMULT)
349        INTEGER BLOCK_CYCLIC_2D , CSRC_ , CTXT_ , DLEN_ , DTYPE_ ,
350       $LLD_ , MB_ , M_ , NB_ , N_ , RSRC_
351        PARAMETER( BLOCK_CYCLIC_2D = 1 , DLEN_ = 9 , DTYPE_ = 1 ,
352       $CTXT_ = 2 , M_ = 3 , N_ = 4 , MB_ = 5 , NB_ = 6 ,
353       $RSRC_ = 7 , CSRC_ = 8 , LLD_ = 9 )
354        CALL CLACPY( 'N' , MAX_BW , MAX_BW , AF( WORK_U + ODD_SIZE*BWL + 1 ) ,
355       $MAX_BW , AF( ODD_SIZE*BWU + 2*MBW2 + 1 ) , MAX_BW )
356  
357        ELSE
358            COMM_PROC = MYCOL + LEVEL_DIST / 2
359        ENDIF
360  
361        IF( MYCOL / LEVEL_DIST .LE.(NPCOL - 1) / LEVEL_DIST - 2 )THEN
362  
363            CALL CGERV2D( ICTXT , MAX_BW , MAX_BW ,
364       $    AF( ODD_SIZE*BWU + 1 ) ,
365       $    MAX_BW , 0 , COMM_PROC )
366  
367            CALL CGERV2D( ICTXT , MAX_BW , MAX_BW ,
368       $    AF( WORK_U + ODD_SIZE*BWL + 1 ) ,
369       $    MAX_BW , 0 , COMM_PROC )
370  
371            IF( INFO .EQ. 0 ) THEN
372  
373  *             Modify upper off_diagonal block with diagonal block
374  
375                CALL CTBTRS(
376       $        'L' , 'N' , 'U' , BWU , MIN( BWL , BWU - 1 ) , BWU ,
377       $        AF( ODD_SIZE*BWU +
378       $        MBW2 + 1 + (MAX_BW + 1)*(MAX_BW - BWU)) , MAX_BW + 1 ,
379       $        AF( WORK_U + ODD_SIZE*BWL + 1 + MAX_BW - BWU) , MAX_BW , INFO )
380  
381  *             Modify lower off_diagonal block with diagonal block
382  
383                CALL CTBTRS(
384       $        'U' , 'C' , 'N' , BWL , MIN( BWU , BWL - 1 ) , BWL ,
385       $        AF( ODD_SIZE*BWU +
386       $        MBW2 + 1 - MIN( BWU , BWL - 1 ) + (MAX_BW + 1)*(MAX_BW - BWL)) , MAX_BW + 1 ,
387       $        AF( ODD_SIZE*BWU + 1 + MAX_BW - BWL) , MAX_BW , INFO )
388  
389            ENDIF
390  *         End of "if( info.eq.0 ) then"
391  
392  *         Calculate contribution from this block to next diagonal block
393  
394            CALL CGEMM( 'C' , 'N' , MAX_BW , MAX_BW , MAX_BW , - CONE ,
395       $    AF((ODD_SIZE)*BWU + 1 ) , MAX_BW ,
396       $    AF( WORK_U + (ODD_SIZE)*BWL + 1 ) , MAX_BW , CZERO ,
397       $    WORK( 1 ) , MAX_BW )
398  
399  *         Send contribution to diagonal block's owning processor.
400  
401            CALL CGESD2D( ICTXT , MAX_BW , MAX_BW , WORK( 1 ) , MAX_BW ,
402       $    0 , MYCOL + LEVEL_DIST )
403  
404        ENDIF
405  *     End of "if( mycol/level_dist .le.(npcol-1)/level_dist-2 )..."
406  
407  *     ****************************************************************
408  *     Receive off_diagonal block from left and use to finish with this
409  *     processor.
410  
411        IF((MYCOL / LEVEL_DIST .GT. 0 ).AND.
412       $( MYCOL / LEVEL_DIST .LE.(NPCOL - 1) / LEVEL_DIST - 1 ) ) THEN
413  
414        IF( LEVEL_DIST .GT. 1)THEN
415  
416  *         Receive offdiagonal block(s) from proc level_dist / 2 to the
417  *         left
418  
419            CALL CGERV2D( ICTXT , MAX_BW , MAX_BW ,
420       $    AF( WORK_U + ODD_SIZE*BWL + 2*MBW2 + 1 ) ,
421       $    MAX_BW , 0 , MYCOL - LEVEL_DIST / 2 )
422  
423  *         Receive offdiagonal block(s) from proc level_dist / 2 to the
424  *         left
425  
426            CALL CGERV2D( ICTXT , MAX_BW , MAX_BW ,
427       $    AF( ODD_SIZE*BWU + 2*MBW2 + 1 ) ,
428       $    MAX_BW , 0 , MYCOL - LEVEL_DIST / 2 )
429  
430        ENDIF
431  
432        IF( INFO.EQ.0 ) THEN
433  
434  *         Use diagonal block(s) to modify this offdiagonal block
435  
436  *         Since CTBTRS has no "left-right" option , we must transpose
437  
438            CALL CLATCPY( 'N' , MAX_BW , MAX_BW , AF( WORK_U + ODD_SIZE*BWL +
439       $    2*MBW2 + 1) , MAX_BW , WORK( 1 ) , MAX_BW )
440  
441            CALL CTBTRS(
442       $    'L' , 'N' , 'U' , MAX_BW , MIN( BWL , MAX_BW - 1 ) , BWL ,
443       $    AF( ODD_SIZE*BWU + MBW2 + 1) , MAX_BW + 1 ,
444       $    WORK( 1 + MAX_BW*(MAX_BW - BWL) ) , MAX_BW , INFO )
445  
446  *         Transpose back
447  
448            CALL CLATCPY(
449       $    'N' , MAX_BW , MAX_BW , WORK( 1 ) , MAX_BW ,
450       $    AF( WORK_U + ODD_SIZE*BWL +
451       $    2*MBW2 + 1) , MAX_BW )
452  
453  *         Since CTBTRS has no "left-right" option , we must transpose
454  
455            CALL CLATCPY( 'N' , MAX_BW , MAX_BW , AF( ODD_SIZE*BWU +
456       $    2*MBW2 + 1) , MAX_BW , WORK( 1 ) , MAX_BW )
457  
458            CALL CTBTRS(
459       $    'U' , 'C' , 'N' , MAX_BW , MIN( BWU , MAX_BW - 1 ) , BWU ,
460       $    AF( ODD_SIZE*BWU + MBW2 + 1 - MIN( BWU , MAX_BW - 1 )) , MAX_BW + 1 ,
461       $    WORK( 1 + MAX_BW*(MAX_BW - BWU) ) , MAX_BW , INFO )
462  
463  *         Transpose back
464  
465            CALL CLATCPY(
466       $    'N' , MAX_BW , MAX_BW , WORK( 1 ) , MAX_BW ,
467       $    AF( ODD_SIZE*BWU +
468       $    2*MBW2 + 1) , MAX_BW )
469  
470        ENDIF
471  *     End of "if( info.eq.0 ) then"
472  
473  *     Use offdiag block(s) to calculate modification to diag block
474  *     of processor to the left
475  
476        CALL CGEMM( 'N' , 'C' , MAX_BW , MAX_BW , MAX_BW , - CONE ,
477       $AF((ODD_SIZE)*BWU + 2*MBW2 + 1 ) , MAX_BW ,
478       $AF( WORK_U + (ODD_SIZE)*BWL + 2*MBW2 + 1 ) , MAX_BW ,
479       $CZERO , WORK( 1 ) , MAX_BW )
480  
481  *     Send contribution to diagonal block's owning processor.
482  
483        CALL CGESD2D( ICTXT , MAX_BW , MAX_BW , WORK( 1 ) , MAX_BW ,
484       $0 , MYCOL - LEVEL_DIST )
485  
486  *     *******************************************************
487  
488        IF( MYCOL / LEVEL_DIST .LE.(NPCOL - 1) / LEVEL_DIST - 2 ) THEN
489  
490  *         Decide which processor offdiagonal block(s) goes to
491  
492            IF(( MOD( MYCOL / ( 2*LEVEL_DIST ) , 2 )) .EQ.0 ) THEN
493                COMM_PROC = MYCOL + LEVEL_DIST
494            ELSE
495                COMM_PROC = MYCOL - LEVEL_DIST
496            ENDIF
497  
498  *         Use offdiagonal blocks to calculate offdiag
499  *         block to send to neighboring processor. Depending
500  *         on circumstances , may need to transpose the matrix.
501  
502            CALL CGEMM( 'N' , 'N' , MAX_BW , MAX_BW , MAX_BW , - CONE ,
503       $    AF( WORK_U + ODD_SIZE*BWL + 2*MBW2 + 1 ) , MAX_BW ,
504       $    AF( ODD_SIZE*BWU + 1 ) , MAX_BW , CZERO ,
505       $    WORK( 1 ) , MAX_BW )
506  
507  *         Send contribution to offdiagonal block's owning processor.
508  
509            CALL CGESD2D( ICTXT , MAX_BW , MAX_BW , WORK( 1 ) , MAX_BW ,
510       $    0 , COMM_PROC )
511  
512            CALL CGEMM( 'N' , 'N' , MAX_BW , MAX_BW , MAX_BW , - CONE ,
513       $    AF( ODD_SIZE*BWU + 2*MBW2 + 1 ) , MAX_BW ,
514       $    AF( WORK_U + ODD_SIZE*BWL + 1 ) , MAX_BW , CZERO ,
515       $    WORK( 1 ) , MAX_BW )
516  
517  *         Send contribution to offdiagonal block's owning processor.
518  
519            CALL CGESD2D( ICTXT , MAX_BW , MAX_BW , WORK( 1 ) , MAX_BW ,
520       $    0 , COMM_PROC )
521  
522        ENDIF
523  
524        ENDIF
525  *     End of "if( mycol/level_dist.le.(npcol-1)/level_dist -1 )..."
526  
527     14 CONTINUE
528  
529   1000 CONTINUE
530  
531  *     Free BLACS space used to hold standard - form grid.
532  
533        IF( ICTXT_SAVE .NE. ICTXT_NEW ) THEN
534            CALL BLACS_GRIDEXIT( ICTXT_NEW )
535        ENDIF
536  
537   1234 CONTINUE
538  
539  *     Restore saved input parameters
540  
541        ICTXT = ICTXT_SAVE
542        NP = NP_SAVE
543  
544  *     Output minimum worksize
545  
546        WORK( 1 ) = WORK_SIZE_MIN
547  
548  *     Make INFO consistent across processors
549  
550        CALL IGAMX2D( ICTXT , 'A' , ' ' , 1 , 1 , INFO , 1 , INFO , INFO ,
551       $- 1 , 0 , 0 )
552  
553        IF( MYCOL.EQ.0 ) THEN
554            CALL IGEBS2D( ICTXT , 'A' , ' ' , 1 , 1 , INFO , 1 )
555        ELSE
556            CALL IGEBR2D( ICTXT , 'A' , ' ' , 1 , 1 , INFO , 1 , 0 , 0 )
557        ENDIF
558  
559        RETURN
560  
561  *     End of PCDBTRF
562  
563        END