Routine: PDPBTRF()  File: SRC\pdpbtrf.f

 
 
# lines: 1475
  # code: 1475
  # comment: 0
  # blank:0
# Variables:35
# Callers:1
# Callings:0
# Words:455
# Keywords:191
 

 

 
..
     .. 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
     Pre-calculate bw^2
     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
       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 {L_i}^T in each processor
         Apply factorization to odd-even connection block B_i
         transpose the connection block in preparation.
         Perform the triangular system solve {L_i}{{B'}_i}^T = {B_i}^T
         transpose resulting block to its location
           in main storage.
         Compute contribution to diagonal block(s) of reduced system.
          {C'}_i = {C_i}-{{B'}_i}{{B'}_i}^T
         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*bw)*bw+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} {{G}_i}^T = {D_i}$ .
         Calculate the update block for previous proc, E_i = G_i{G_i}^T
         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 transpose form.
           {F_i}^T =  {H_i}{{B'}_i}^T
           Copy matrix H_i (the last bw cols of G_i) to AF storage
             as per requirements of BLAS routine DTRMM.
             Since we have G_i^T stored, transpose
             H_i^T to H_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.
         Modify upper off_diagonal block with diagonal block
         End of "if ( info.eq.0 ) then"
         Calculate contribution from this block to next diagonal block
         Send contribution to diagonal block's owning processor.
       End of "if( mycol/level_dist .le. (npcol-1)/level_dist-2 )..."
       ****************************************************************
       Receive off_diagonal block from left and use to finish with this
         processor.
           Receive offdiagonal block(s) from proc level_dist/2 to the
           left
 
Display dynamic version Find AutoScroll Reload FontSize: - + Hide Comments Hide Blanks Frame FullScreen MailPrint

  
001        SUBROUTINE PDPBTRF( UPLO , N , BW , 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  *     April 3 , 2000
008  
009  *     .. Scalar Arguments ..
010        CHARACTER UPLO
011        INTEGER BW , INFO , JA , LAF , LWORK , N
012  *     ..
013  *     .. Array Arguments ..
014        INTEGER DESCA( * )
015        DOUBLE PRECISION A( * ) , AF( * ) , WORK( * )
016  *     ..
017  
018  *     Purpose
019  *     === ====
020  
021  *     PDPBTRF computes a Cholesky factorization
022  *     of an N - by - N real banded
023  *     symmetric positive definite distributed matrix
024  *     with bandwidth BW : A(1 : N , JA : JA + N - 1).
025  *     Reordering is used to increase parallelism in the factorization.
026  *     This reordering results in factors that are DIFFERENT from those
027  *     produced by equivalent sequential codes. These factors cannot
028  *     be used directly by users ; however , they can be used in
029  *     subsequent calls to PDPBTRS to solve linear systems.
030  
031  *     The factorization has the form
032  
033  *     P A(1 : N , JA : JA + N - 1) P^T = U' U , if UPLO = 'U' , or
034  
035  *     P A(1 : N , JA : JA + N - 1) P^T = L L' , if UPLO = 'L'
036  
037  *     where U is a banded upper triangular matrix and L is banded
038  *     lower triangular , and P is a permutation matrix.
039  
040  *     === ==================================================================
041  
042  *     Arguments
043  *     === ======
044  
045  *     UPLO(global input) CHARACTER
046  *     = 'U' : Upper triangle of A(1 : N , JA : JA + N - 1) is stored ;
047  *     = 'L' : Lower triangle of A(1 : N , JA : JA + N - 1) is stored.
048  
049  *     N(global input) INTEGER
050  *     The number of rows and columns to be operated on , i.e. the
051  *     order of the distributed submatrix A(1 : N , JA : JA + N - 1). N >= 0.
052  
053  *     BW(global input) INTEGER
054  *     Number of subdiagonals in L or U. 0 <= BW <= N - 1
055  
056  *     A(local input / local output) DOUBLE PRECISION pointer into
057  *     local memory to an array with first dimension
058  *     LLD_A >=(bw + 1)(stored in DESCA).
059  *     On entry , this array contains the local pieces of the
060  *     N - by - N symmetric banded distributed matrix
061  *     A(1 : N , JA : JA + N - 1) to be factored.
062  *     This local portion is stored in the packed banded format
063  *     used in LAPACK. Please see the Notes below and the
064  *     ScaLAPACK manual for more detail on the format of
065  *     distributed matrices.
066  *     On exit , this array contains information containing details
067  *     of the factorization.
068  *     Note that permutations are performed on the matrix , so that
069  *     the factors returned are different from those returned
070  *     by LAPACK.
071  
072  *     JA(global input) INTEGER
073  *     The index in the global array A that points to the start of
074  *     the matrix to be operated on(which may be either all of A
075  *     or a submatrix of A).
076  
077  *     DESCA(global and local input) INTEGER array of dimension DLEN.
078  *     if 1D type(DTYPE_A = 501) , DLEN >= 7 ;
079  *     if 2D type(DTYPE_A = 1) , DLEN >= 9 .
080  *     The array descriptor for the distributed matrix A.
081  *     Contains information of mapping of A to memory. Please
082  *     see NOTES below for full description and options.
083  
084  *     AF(local output) DOUBLE PRECISION array , dimension LAF.
085  *     Auxiliary Fillin Space.
086  *     Fillin is created during the factorization routine
087  *     PDPBTRF and this is stored in AF. If a linear system
088  *     is to be solved using PDPBTRS after the factorization
089  *     routine , AF *must not be altered* after the factorization.
090  
091  *     LAF(local input) INTEGER
092  *     Size of user - input Auxiliary Fillin space AF. Must be >=
093  *     (NB + 2*bw)*bw
094  *     If LAF is not large enough , an error code will be returned
095  *     and the minimum acceptable size will be returned in AF( 1 )
096  
097  *     WORK(local workspace / local output)
098  *     DOUBLE PRECISION temporary workspace. This space may
099  *     be overwritten in between calls to routines. WORK must be
100  *     the size given in LWORK.
101  *     On exit , WORK( 1 ) contains the minimal LWORK.
102  
103  *     LWORK(local input or global input) INTEGER
104  *     Size of user - input workspace WORK.
105  *     If LWORK is too small , the minimal acceptable size will be
106  *     returned in WORK(1) and an error code is returned. LWORK >=
107  *     bw*bw
108  
109  *     INFO(global output) INTEGER
110  *     = 0 : successful exit
111  *     < 0 : If the i - th argument is an array and the j - entry had
112  *     an illegal value , then INFO = - (i*100 + j) , if the i - th
113  *     argument is a scalar and had an illegal value , then
114  *     INFO = - i.
115  *     > 0 : If INFO = K <= NPROCS , the submatrix stored on processor
116  *     INFO and factored locally was not
117  *     positive definite , and
118  *     the factorization was not completed.
119  *     If INFO = K > NPROCS , the submatrix stored on processor
120  *     INFO - NPROCS representing interactions with other
121  *     processors was not
122  *     positive definite ,
123  *     and the factorization was not completed.
124  
125  *     === ==================================================================
126  
127  *     Restrictions
128  *     === =========
129  
130  *     The following are restrictions on the input parameters. Some of these
131  *     are temporary and will be removed in future releases , while others
132  *     may reflect fundamental technical limitations.
133  
134  *     Non - cyclic restriction : VERY IMPORTANT !
135  *     P*NB >= mod(JA - 1 , NB) + N.
136  *     The mapping for matrices must be blocked , reflecting the nature
137  *     of the divide and conquer algorithm as a task - parallel algorithm.
138  *     This formula in words is : no processor may have more than one
139  *     chunk of the matrix.
140  
141  *     Blocksize cannot be too small :
142  *     If the matrix spans more than one processor , the following
143  *     restriction on NB , the size of each block on each processor ,
144  *     must hold :
145  *     NB >= 2*BW
146  *     The bulk of parallel computation is done on the matrix of size
147  *     O(NB) on each processor. If this is too small , divide and conquer
148  *     is a poor choice of algorithm.
149  
150  *     Submatrix reference :
151  *     JA = IB
152  *     Alignment restriction that prevents unnecessary communication.
153  
154  *     === ==================================================================
155  
156  *     Notes
157  *     === ==
158  
159  *     If the factorization routine and the solve routine are to be called
160  *     separately(to solve various sets of righthand sides using the same
161  *     coefficient matrix) , the auxiliary space AF *must not be altered*
162  *     between calls to the factorization routine and the solve routine.
163  
164  *     The best algorithm for solving banded and tridiagonal linear systems
165  *     depends on a variety of parameters , especially the bandwidth.
166  *     Currently , only algorithms designed for the case N / P >> bw are
167  *     implemented. These go by many names , including Divide and Conquer ,
168  *     Partitioning , domain decomposition - type , etc.
169  
170  *     Algorithm description : Divide and Conquer
171  
172  *     The Divide and Conqer algorithm assumes the matrix is narrowly
173  *     banded compared with the number of equations. In this situation ,
174  *     it is best to distribute the input matrix A one - dimensionally ,
175  *     with columns atomic and rows divided amongst the processes.
176  *     The basic algorithm divides the banded matrix up into
177  *     P pieces with one stored on each processor ,
178  *     and then proceeds in 2 phases for the factorization or 3 for the
179  *     solution of a linear system.
180  *     1) Local Phase :
181  *     The individual pieces are factored independently and in
182  *     parallel. These factors are applied to the matrix creating
183  *     fillin , which is stored in a non - inspectable way in auxiliary
184  *     space AF. Mathematically , this is equivalent to reordering
185  *     the matrix A as P A P^T and then factoring the principal
186  *     leading submatrix of size equal to the sum of the sizes of
187  *     the matrices factored on each processor. The factors of
188  *     these submatrices overwrite the corresponding parts of A
189  *     in memory.
190  *     2) Reduced System Phase :
191  *     A small(BW* (P - 1)) system is formed representing
192  *     interaction of the larger blocks , and is stored(as are its
193  *     factors) in the space AF. A parallel Block Cyclic Reduction
194  *     algorithm is used. For a linear system , a parallel front solve
195  *     followed by an analagous backsolve , both using the structure
196  *     of the factored matrix , are performed.
197  *     3) Backsubsitution Phase :
198  *     For a linear system , a local backsubstitution is performed on
199  *     each processor in parallel.
200  
201  *     Descriptors
202  *     === ========
203  
204  *     Descriptors now have *types* and differ from ScaLAPACK 1.0.
205  
206  *     Note : banded codes can use either the old two dimensional
207  *     or new one - dimensional descriptors , though the processor grid in
208  *     both cases *must be one - dimensional*. We describe both types below.
209  
210  *     Each global data object is described by an associated description
211  *     vector. This vector stores the information required to establish
212  *     the mapping between an object element and its corresponding process
213  *     and memory location.
214  
215  *     Let A be a generic term for any 2D block cyclicly distributed array.
216  *     Such a global array has an associated description vector DESCA.
217  *     In the following comments , the character _ should be read as
218  *     "of the global array".
219  
220  *     NOTATION STORED IN EXPLANATION
221  *     --- ------------ -------------- --------------------------------------
222  *     DTYPE_A(global) DESCA( DTYPE_ )The descriptor type. In this case ,
223  *     DTYPE_A = 1.
224  *     CTXT_A(global) DESCA( CTXT_ ) The BLACS context handle , indicating
225  *     the BLACS process grid A is distribu -
226  *     ted over. The context itself is glo -
227  *     bal , but the handle(the integer
228  *     value) may vary.
229  *     M_A(global) DESCA( M_ ) The number of rows in the global
230  *     array A.
231  *     N_A(global) DESCA( N_ ) The number of columns in the global
232  *     array A.
233  *     MB_A(global) DESCA( MB_ ) The blocking factor used to distribute
234  *     the rows of the array.
235  *     NB_A(global) DESCA( NB_ ) The blocking factor used to distribute
236  *     the columns of the array.
237  *     RSRC_A(global) DESCA( RSRC_ ) The process row over which the first
238  *     row of the array A is distributed.
239  *     CSRC_A(global) DESCA( CSRC_ ) The process column over which the
240  *     first column of the array A is
241  *     distributed.
242  *     LLD_A(local) DESCA( LLD_ ) The leading dimension of the local
243  *     array. LLD_A >= MAX(1 , LOCr(M_A)).
244  
245  *     Let K be the number of rows or columns of a distributed matrix ,
246  *     and assume that its process grid has dimension p x q.
247  *     LOCr( K ) denotes the number of elements of K that a process
248  *     would receive if K were distributed over the p processes of its
249  *     process column.
250  *     Similarly , LOCc( K ) denotes the number of elements of K that a
251  *     process would receive if K were distributed over the q processes of
252  *     its process row.
253  *     The values of LOCr() and LOCc() may be determined via a call to the
254  *     ScaLAPACK tool function , NUMROC :
255  *     LOCr( M ) = NUMROC( M , MB_A , MYROW , RSRC_A , NPROW ) ,
256  *     LOCc( N ) = NUMROC( N , NB_A , MYCOL , CSRC_A , NPCOL ).
257  *     An upper bound for these quantities may be computed by :
258  *     LOCr( M ) <= ceil( ceil(M / MB_A) / NPROW )*MB_A
259  *     LOCc( N ) <= ceil( ceil(N / NB_A) / NPCOL )*NB_A
260  
261  *     One - dimensional descriptors :
262  
263  *     One - dimensional descriptors are a new addition to ScaLAPACK since
264  *     version 1.0. They simplify and shorten the descriptor for 1D
265  *     arrays.
266  
267  *     Since ScaLAPACK supports two - dimensional arrays as the fundamental
268  *     object , we allow 1D arrays to be distributed either over the
269  *     first dimension of the array(as if the grid were P - by - 1) or the
270  *     2nd dimension(as if the grid were 1 - by - P). This choice is
271  *     indicated by the descriptor type(501 or 502)
272  *     as described below.
273  
274  *     IMPORTANT NOTE : the actual BLACS grid represented by the
275  *     CTXT entry in the descriptor may be *either* P - by - 1 or 1 - by - P
276  *     irrespective of which one - dimensional descriptor type
277  *     (501 or 502) is input.
278  *     This routine will interpret the grid properly either way.
279  *     ScaLAPACK routines *do not support intercontext operations* so that
280  *     the grid passed to a single ScaLAPACK routine *must be the same*
281  *     for all array descriptors passed to that routine.
282  
283  *     NOTE : In all cases where 1D descriptors are used , 2D descriptors
284  *     may also be used , since a one - dimensional array is a special case
285  *     of a two - dimensional array with one dimension of size unity.
286  *     The two - dimensional array used in this case *must* be of the
287  *     proper orientation :
288  *     If the appropriate one - dimensional descriptor is DTYPEA = 501
289  *     (1 by P type) , then the two dimensional descriptor must
290  *     have a CTXT value that refers to a 1 by P BLACS grid ;
291  *     If the appropriate one - dimensional descriptor is DTYPEA = 502
292  *     (P by 1 type) , then the two dimensional descriptor must
293  *     have a CTXT value that refers to a P by 1 BLACS grid.
294  
295  *     Summary of allowed descriptors , types , and BLACS grids :
296  *     DTYPE 501 502 1 1
297  *     BLACS grid 1xP or Px1 1xP or Px1 1xP Px1
298  *     --- --------------------------------------------------
299  *     A               OK NO OK NO
300  *     B               NO OK NO OK
301  
302  *     Let A be a generic term for any 1D block cyclicly distributed array.
303  *     Such a global array has an associated description vector DESCA.
304  *     In the following comments , the character _ should be read as
305  *     "of the global array".
306  
307  *     NOTATION STORED IN EXPLANATION
308  *     --- ------------ ---------- ------------------------------------------
309  *     DTYPE_A(global) DESCA( 1 ) The descriptor type. For 1D grids ,
310  *     TYPE_A = 501 : 1 - by - P grid.
311  *     TYPE_A = 502 : P - by - 1 grid.
312  *     CTXT_A(global) DESCA( 2 ) The BLACS context handle , indicating
313  *     the BLACS process grid A is distribu -
314  *     ted over. The context itself is glo -
315  *     bal , but the handle(the integer
316  *     value) may vary.
317  *     N_A(global) DESCA( 3 ) The size of the array dimension being
318  *     distributed.
319  *     NB_A(global) DESCA( 4 ) The blocking factor used to distribute
320  *     the distributed dimension of the array.
321  *     SRC_A(global) DESCA( 5 ) The process row or column over which the
322  *     first row or column of the array
323  *     is distributed.
324  *     LLD_A(local) DESCA( 6 ) The leading dimension of the local array
325  *     storing the local blocks of the distri -
326  *     buted array A. Minimum value of LLD_A
327  *     depends on TYPE_A.
328  *     TYPE_A = 501 : LLD_A >=
329  *     size of undistributed dimension , 1.
330  *     TYPE_A = 502 : LLD_A >= NB_A , 1.
331  *     Reserved DESCA( 7 ) Reserved for future use.
332  
333  *     === ==================================================================
334  
335  *     Code Developer : Andrew J. Cleary , University of Tennessee.
336  *     Current address : Lawrence Livermore National Labs.
337  
338  *     === ==================================================================
339  
340  *     .. Parameters ..
341        DOUBLE PRECISION ONE
342        PARAMETER( ONE = 1.0D + 0 )
343        DOUBLE PRECISION ZERO
344        PARAMETER( ZERO = 0.0D + 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        END IF
355  
356        IF( INFO.EQ.0 ) THEN
357  
358  *         Use diagonal block(s) to modify this offdiagonal block
359  
360            CALL DTRSM( 'R' , 'L' , 'T' , 'N' , BW , BW , ONE ,
361       $    AF( ODD_SIZE*BW + MBW2 + 1 ) , BW ,
362       $    AF( ODD_SIZE*BW + 2*MBW2 + 1 ) , BW )
363  
364        END IF
365  *     End of "if( info.eq.0 ) then"
366  
367  *     Use offdiag block(s) to calculate modification to diag block
368  *     of processor to the left
369  
370        CALL DSYRK( 'L' , 'N' , BW , BW , - ONE ,
371       $AF(( ODD_SIZE + 2*BW )*BW + 1 ) , BW , ZERO ,
372       $WORK( 1 ) , BW )
373  
374  *     Send contribution to diagonal block's owning processor.
375  
376        CALL DGESD2D( ICTXT , BW , BW , WORK( 1 ) , BW , 0 ,
377       $MYCOL - LEVEL_DIST )
378  
379  *     *******************************************************
380  
381        IF( MYCOL / LEVEL_DIST.LE.( NPCOL - 1 ) / LEVEL_DIST - 2 ) THEN
382  
383  *         Decide which processor offdiagonal block(s) goes to
384  
385            IF(( MOD( MYCOL / ( 2*LEVEL_DIST ) , 2 ) ).EQ.0 ) THEN
386                COMM_PROC = MYCOL + LEVEL_DIST
387            ELSE
388                COMM_PROC = MYCOL - LEVEL_DIST
389            END IF
390  
391  *         Use offdiagonal blocks to calculate offdiag
392  *         block to send to neighboring processor. Depending
393  *         on circumstances , may need to transpose the matrix.
394  
395            CALL DGEMM( 'N' , 'N' , BW , BW , BW , - ONE ,
396       $    AF( ODD_SIZE*BW + 2*MBW2 + 1 ) , BW ,
397       $    AF( ODD_SIZE*BW + 1 ) , BW , ZERO , WORK( 1 ) ,
398       $    BW )
399  
400  *         Send contribution to offdiagonal block's owning processor.
401  
402            CALL DGESD2D( ICTXT , BW , BW , WORK( 1 ) , BW , 0 ,
403       $    COMM_PROC )
404  
405        END IF
406  
407        END IF
408  *     End of "if( mycol/level_dist.le.(npcol-1)/level_dist -1 )..."
409  
410     60 CONTINUE
411  
412        ELSE
413  
414  *         CASE UPLO = 'U'
415  
416  *         *******************************************************************
417  *         PHASE 1 : Local computation phase.
418  *         *******************************************************************
419  
420  *         Sizes of the extra triangles communicated bewtween processors
421  
422            IF( MYCOL.GT.0 ) THEN
423                PREV_TRI_SIZE_M = MIN( BW , NUMROC( N , PART_SIZE , MYCOL , 0 ,
424       $        NPCOL ) )
425                PREV_TRI_SIZE_N = MIN( BW , NUMROC( N , PART_SIZE , MYCOL - 1 , 0 ,
426       $        NPCOL ) )
427            END IF
428  
429            IF( MYCOL.LT.NPCOL - 1 ) THEN
430                NEXT_TRI_SIZE_M = MIN( BW , NUMROC( N , PART_SIZE , MYCOL + 1 , 0 ,
431       $        NPCOL ) )
432                NEXT_TRI_SIZE_N = MIN( BW , NUMROC( N , PART_SIZE , MYCOL , 0 ,
433       $        NPCOL ) )
434            END IF
435  
436  *         Factor main partition A_i^T = U_i {U_i}^T in each processor
437  
438            CALL DPBTRF( UPLO , ODD_SIZE , BW , A( OFST + 1 ) , LLDA , INFO )
439  
440            IF( INFO.NE.0 ) THEN
441                INFO = MYCOL + 1
442                GO TO 70
443            END IF
444  
445            IF( MYCOL.LT.NP - 1 ) THEN
446  *             Apply factorization to odd - even connection block B_i
447  
448  *             Move the connection block in preparation.
449  
450                CALL DLACPY( 'L' , BW , BW , A(( OFST + 1 + ODD_SIZE*LLDA ) ) ,
451       $        LLDA - 1 , AF( ODD_SIZE*BW + 2*MBW2 + 1 + BW - BW ) , BW )
452  
453  *             Perform the triangular solve {L_i}{{B'}_i}^T = {B_i}^T
454  
455                CALL DTRTRS( 'U' , 'T' , 'N' , BW , BW ,
456       $        A( OFST + BW + 1 + ( ODD_SIZE - BW )*LLDA ) , LLDA - 1 ,
457       $        AF( ODD_SIZE*BW + 2*MBW2 + 1 ) , BW , INFO )
458  
459  *             Move the resulting block back to its location in main storage.
460  
461                CALL DLACPY( 'L' , BW , BW , AF( ODD_SIZE*BW + 2*MBW2 + 1 + BW - BW ) ,
462       $        BW , A(( OFST + 1 + ODD_SIZE*LLDA ) ) , LLDA - 1 )
463  
464  *             Compute contribution to diagonal block(s) of reduced system.
465  *             {C'}_i^T = {C_i}^T - {{B'}_i}^T{{B'}_i}
466  
467  *             The following method uses more flops than necessary but
468  *             does not necessitate the writing of a new BLAS routine.
469  
470                CALL DSYRK( UPLO , 'T' , BW , BW , - ONE ,
471       $        AF( ODD_SIZE*BW + 2*MBW2 + 1 ) , BW , ONE ,
472       $        A( OFST + BW + 1 + ODD_SIZE*LLDA ) , LLDA - 1 )
473  
474            END IF
475  *         End of "if( MYCOL .lt. NP-1 )..." loop
476  
477     70 CONTINUE
478  *     If the processor could not locally factor , it jumps here.
479  
480        IF( MYCOL.NE.0 ) THEN
481  *         Discard temporary matrix stored beginning in
482  *         AF((odd_size + 2*bw)*bw + 1 ) and use for
483  *         off_diagonal block of reduced system.
484  
485  *         Calculate the "spike" fillin , ${L_i} {{G}_i}^T = {D_i}$ .
486  
487  *         Copy D block into AF storage for solve.
488  
489            CALL DLATCPY( 'L' , PREV_TRI_SIZE_N , PREV_TRI_SIZE_M ,
490       $    A( OFST + 1 ) , LLDA - 1 , AF( 1 ) , ODD_SIZE )
491  
492            IF( INFO.EQ.0 ) THEN
493  
494                CALL DTBTRS( 'U' , 'T' , 'N' , ODD_SIZE , BW , BW ,
495       $        A( OFST + 1 ) , LLDA , AF( 1 ) , ODD_SIZE , INFO )
496  
497  *             Calculate the update block for previous proc , E_i = G_i{G_i}^T
498  
499                CALL DSYRK( 'L' , 'T' , BW , ODD_SIZE , - ONE , AF( 1 ) ,
500       $        ODD_SIZE , ZERO , AF( 1 + ( ODD_SIZE + 2*BW )*BW ) ,
501       $        BW )
502  
503  *             Initiate send of E_i to previous processor to overlap
504  *             with next computation.
505  
506                CALL DGESD2D( ICTXT , BW , BW , AF( ODD_SIZE*BW + 2*MBW2 + 1 ) ,
507       $        BW , 0 , MYCOL - 1 )
508  
509                IF( MYCOL.LT.NP - 1 ) THEN
510  
511  *                 Calculate off - diagonal block(s) of reduced system.
512  *                 Note : for ease of use in solution of reduced system , store
513  *                 L's off - diagonal block in transpose form.
514  *                 {F_i}^T = {H_i}{{B'}_i}^T
515  
516  *                 Copy matrix H_i(the last bw cols of G_i) to AF storage
517  *                 as per requirements of BLAS routine DTRMM.
518  *                 Since we have G_i^T stored , transpose
519  *                 H_i^T to H_i.
520  
521                    CALL DLATCPY( 'N' , BW , BW , AF( ODD_SIZE - BW + 1 ) ,
522       $            ODD_SIZE , AF(( ODD_SIZE )*BW + 1 ) , BW )
523  
524                    CALL DTRMM( 'R' , 'L' , 'N' , 'N' , BW , BW , - ONE ,
525       $            A(( OFST + 1 + ODD_SIZE*LLDA ) ) , LLDA - 1 ,
526       $            AF(( ODD_SIZE )*BW + 1 ) , BW )
527  
528                END IF
529  
530            END IF
531  *         End of "if( MYCOL .ne. 0 )..."
532  
533        END IF
534  *     End of "if(info.eq.0) then"
535  
536  *     Check to make sure no processors have found errors
537  
538        CALL IGAMX2D( ICTXT , 'A' , ' ' , 1 , 1 , INFO , 1 , INFO , INFO , - 1 ,
539       $0 , 0 )
540  
541        IF( MYCOL.EQ.0 ) THEN
542            CALL IGEBS2D( ICTXT , 'A' , ' ' , 1 , 1 , INFO , 1 )
543        ELSE
544            CALL IGEBR2D( ICTXT , 'A' , ' ' , 1 , 1 , INFO , 1 , 0 , 0 )
545        END IF
546  
547        IF( INFO.NE.0 ) THEN
548            GO TO 110
549        END IF
550  *     No errors found , continue
551  
552  *     *******************************************************************
553  *     PHASE 2 : Formation and factorization of Reduced System.
554  *     *******************************************************************
555  
556  *     Gather up local sections of reduced system
557  
558  *     The last processor does not participate in the factorization of
559  *     the reduced system , having sent its E_i already.
560        IF( MYCOL.EQ.NPCOL - 1 ) THEN
561            GO TO 100
562        END IF
563  
564  *     Initiate send of off - diag block(s) to overlap with next part.
565  *     Off - diagonal block needed on neighboring processor to start
566  *     algorithm.
567  
568        IF(( MOD( MYCOL + 1 , 2 ).EQ.0 ) .AND.( MYCOL.GT.0 ) ) THEN
569  
570            CALL DGESD2D( ICTXT , BW , BW , AF( ODD_SIZE*BW + 1 ) , BW , 0 ,
571       $    MYCOL - 1 )
572  
573        END IF
574  
575  *     Transpose last diagonal block into AF storage for subsequent
576  *     operations.
577  
578        CALL DLATCPY( 'U' , BW , BW , A( OFST + ODD_SIZE*LLDA + 1 + BW ) ,
579       $LLDA - 1 , AF( ODD_SIZE*BW + MBW2 + 1 ) , BW )
580  
581  *     Receive cont. to diagonal block that is stored on this proc.
582  
583        IF( MYCOL.LT.NPCOL - 1 ) THEN
584  
585            CALL DGERV2D( ICTXT , BW , BW , AF( ODD_SIZE*BW + 2*MBW2 + 1 ) , BW ,
586       $    0 , MYCOL + 1 )
587  
588  *         Add contribution to diagonal block
589  
590            CALL DAXPY( MBW2 , ONE , AF( ODD_SIZE*BW + 2*MBW2 + 1 ) , 1 ,
591       $    AF( ODD_SIZE*BW + MBW2 + 1 ) , 1 )
592  
593        END IF
594  
595  *     *************************************
596  *     Modification Loop
597  
598  *     The distance for sending and receiving for each level starts
599  *     at 1 for the first level.
600        LEVEL_DIST = 1
601  
602  *     Do until this proc is needed to modify other procs' equations
603  
604     80 CONTINUE
605        IF( MOD(( MYCOL + 1 ) / LEVEL_DIST , 2 ).NE.0 )
606       $    GO TO 90
607  
608  *         Receive and add contribution to diagonal block from the left
609  
610            IF( MYCOL - LEVEL_DIST.GE.0 ) THEN
611                CALL DGERV2D( ICTXT , BW , BW , WORK( 1 ) , BW , 0 ,
612       $        MYCOL - LEVEL_DIST )
613  
614                CALL DAXPY( MBW2 , ONE , WORK( 1 ) , 1 ,
615       $        AF( ODD_SIZE*BW + MBW2 + 1 ) , 1 )
616  
617            END IF
618  
619  *         Receive and add contribution to diagonal block from the right
620  
621            IF( MYCOL + LEVEL_DIST.LT.NPCOL - 1 ) THEN
622                CALL DGERV2D( ICTXT , BW , BW , WORK( 1 ) , BW , 0 ,
623       $        MYCOL + LEVEL_DIST )
624  
625                CALL DAXPY( MBW2 , ONE , WORK( 1 ) , 1 ,
626       $        AF( ODD_SIZE*BW + MBW2 + 1 ) , 1 )
627  
628            END IF
629  
630            LEVEL_DIST = LEVEL_DIST*2
631  
632            GO TO 80
633     90 CONTINUE
634  *     [End of GOTO Loop]
635  
636  *     *********************************
637  *     Calculate and use this proc's blocks to modify other procs'...
638  
639  *     Factor diagonal block
640  
641        CALL DPOTRF( 'L' , BW , AF( ODD_SIZE*BW + MBW2 + 1 ) , BW , INFO )
642  
643        IF( INFO.NE.0 ) THEN
644            INFO = NPCOL + MYCOL
645        END IF
646  
647  *     ****************************************************************
648  *     Receive offdiagonal block from processor to right.
649  *     If this is the first group of processors , the receive comes
650  *     from a different processor than otherwise.
651  
652        IF( LEVEL_DIST.EQ.1 ) THEN
653            COMM_PROC = MYCOL + 1
654  
655  *         Move block into place that it will be expected to be for
656  *         calcs.
657  
658            CALL DLACPY( 'N' , BW , BW , AF( ODD_SIZE*BW + 1 ) , BW ,
659       $    AF( ODD_SIZE*BW + 2*MBW2 + 1 ) , BW )
660  
661        ELSE
662            COMM_PROC = MYCOL + LEVEL_DIST / 2
663        END IF
664  
665        IF( MYCOL / LEVEL_DIST.LE.( NPCOL - 1 ) / LEVEL_DIST - 2 ) THEN
666  
667            CALL DGERV2D( ICTXT , BW , BW , AF( ODD_SIZE*BW + 1 ) , BW , 0 ,
668       $    COMM_PROC )
669  
670            IF( INFO.EQ.0 ) THEN
671  
672  *             Modify upper off_diagonal block with diagonal block
673  
674                CALL DTRSM( 'L' , 'L' , 'N' , 'N' , BW , BW , ONE ,
675       $        AF( ODD_SIZE*BW + MBW2 + 1 ) , BW ,
676       $        AF( ODD_SIZE*BW + 1 ) , BW )
677  
678            END IF
679  *         End of "if( info.eq.0 ) then"
680  
681  *         Calculate contribution from this block to next diagonal block
682  
683            CALL DSYRK( 'L' , 'T' , BW , BW , - ONE , AF(( ODD_SIZE )*BW + 1 ) ,
684       $    BW , ZERO , WORK( 1 ) , BW )
685  
686  *         Send contribution to diagonal block's owning processor.
687  
688            CALL DGESD2D( ICTXT , BW , BW , WORK( 1 ) , BW , 0 ,
689       $    MYCOL + LEVEL_DIST )
690  
691        END IF
692  *     End of "if( mycol/level_dist .le.(npcol-1)/level_dist-2 )..."
693  
694  *     ****************************************************************
695  *     Receive off_diagonal block from left and use to finish with this
696  *     processor.
697  
698        IF(( MYCOL / LEVEL_DIST.GT.0 ) .AND.
699       $( MYCOL / LEVEL_DIST.LE.( NPCOL - 1 ) / LEVEL_DIST - 1 ) ) THEN
700  
701        IF( LEVEL_DIST.GT.1 ) THEN
702  
703  *         Receive offdiagonal block(s) from proc level_dist / 2 to the
704  *         left
705  
706            CALL DGERV2D( ICTXT , BW , BW , AF( ODD_SIZE*BW + 2*MBW2 + 1 ) ,
707       $    BW , 0 , MYCOL - LEVEL_DIST / 2 )
708  
709        END IF
710  
711        IF( INFO.EQ.0 ) THEN
712  
713  *         Use diagonal block(s) to modify this offdiagonal block
714  
715            CALL DTRSM( 'R' , 'L' , 'T' , 'N' , BW , BW , ONE ,
716       $    AF( ODD_SIZE*BW + MBW2 + 1 ) , BW ,
717       $    AF( ODD_SIZE*BW + 2*MBW2 + 1 ) , BW )
718  
719        END IF
720  *     End of "if( info.eq.0 ) then"
721  
722  *     Use offdiag block(s) to calculate modification to diag block
723  *     of processor to the left
724  
725        CALL DSYRK( 'L' , 'N' , BW , BW , - ONE ,
726       $AF(( ODD_SIZE + 2*BW )*BW + 1 ) , BW , ZERO ,
727       $WORK( 1 ) , BW )
728  
729  *     Send contribution to diagonal block's owning processor.
730  
731        CALL DGESD2D( ICTXT , BW , BW , WORK( 1 ) , BW , 0 ,
732       $MYCOL - LEVEL_DIST )
733  
734  *     *******************************************************
735  
736        IF( MYCOL / LEVEL_DIST.LE.( NPCOL - 1 ) / LEVEL_DIST - 2 ) THEN
737  
738  *         Decide which processor offdiagonal block(s) goes to
739  
740            IF(( MOD( MYCOL / ( 2*LEVEL_DIST ) , 2 ) ).EQ.0 ) THEN
741                COMM_PROC = MYCOL + LEVEL_DIST
742            ELSE
743                COMM_PROC = MYCOL - LEVEL_DIST
744            END IF
745  
746  *         Use offdiagonal blocks to calculate offdiag
747  *         block to send to neighboring processor. Depending
748  *         on circumstances , may need to transpose the matrix.
749  
750            CALL DGEMM( 'N' , 'N' , BW , BW , BW , - ONE ,
751       $    AF( ODD_SIZE*BW + 2*MBW2 + 1 ) , BW ,
752       $    AF( ODD_SIZE*BW + 1 ) , BW , ZERO , WORK( 1 ) ,
753       $    BW )
754  
755  *         Send contribution to offdiagonal block's owning processor.
756  
757            CALL DGESD2D( ICTXT , BW , BW , WORK( 1 ) , BW , 0 ,
758       $    COMM_PROC )
759  
760        END IF
761  
762        END IF
763  *     End of "if( mycol/level_dist.le.(npcol-1)/level_dist -1 )..."
764  
765    100 CONTINUE
766  
767        END IF
768  
769    110 CONTINUE
770  
771  *     Free BLACS space used to hold standard - form grid.
772  
773        IF( ICTXT_SAVE.NE.ICTXT_NEW ) THEN
774            CALL BLACS_GRIDEXIT( ICTXT_NEW )
775        END IF
776  
777    120 CONTINUE
778  
779  *     Restore saved input parameters
780  
781        ICTXT = ICTXT_SAVE
782        NP = NP_SAVE
783  
784  *     Output minimum worksize
785  
786        WORK( 1 ) = WORK_SIZE_MIN
787  
788  *     Make INFO consistent across processors
789  
790        CALL IGAMX2D( ICTXT , 'A' , ' ' , 1 , 1 , INFO , 1 , INFO , INFO , - 1 , 0 ,
791       $0 )
792  
793        IF( MYCOL.EQ.0 ) THEN
794            CALL IGEBS2D( ICTXT , 'A' , ' ' , 1 , 1 , INFO , 1 )
795        ELSE
796            CALL IGEBR2D( ICTXT , 'A' , ' ' , 1 , 1 , INFO , 1 , 0 , 0 )
797        END IF
798  
799        RETURN
800  
801  *     End of PDPBTRF
802  
803        END