Routine: PCPBTRF()  File: SRC\pcpbtrf.f

 
 
# lines: 1510
  # code: 1510
  # comment: 0
  # blank:0
# Variables:37
# Callers:1
# Callings:0
# Words:417
# Keywords:154
 

 

 
..
     .. 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}^C in each processor
         Apply factorization to odd-even connection block B_i
         conjugate transpose the connection block in preparation.
         Perform the triangular system solve {L_i}{{B'}_i}^C = {B_i}^C
         conjugate 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}^C
         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}^C = {D_i}$ .
         Calculate the update block for previous proc, E_i = G_i{G_i}^C
         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.
           {F_i}^C =  {H_i}{{B'}_i}^C
           Copy matrix H_i (the last bw cols of G_i) to AF storage
             as per requirements of BLAS routine CTRMM.
             Since we have G_i^C stored, conjugate transpose
             H_i^C 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 PCPBTRF( 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  *     May 25 , 2001
008  
009  *     .. Scalar Arguments ..
010        CHARACTER UPLO
011        INTEGER BW , INFO , JA , LAF , LWORK , N
012  *     ..
013  *     .. Array Arguments ..
014        INTEGER DESCA( * )
015        COMPLEX A( * ) , AF( * ) , WORK( * )
016  *     ..
017  
018  *     Purpose
019  *     === ====
020  
021  *     PCPBTRF computes a Cholesky factorization
022  *     of an N - by - N complex 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 PCPBTRS 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) COMPLEX 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) COMPLEX array , dimension LAF.
085  *     Auxiliary Fillin Space.
086  *     Fillin is created during the factorization routine
087  *     PCPBTRF and this is stored in AF. If a linear system
088  *     is to be solved using PCPBTRS 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  *     COMPLEX 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  *     This version released : August , 2001.
338  
339  *     === ==================================================================
340  
341  *     ..
342  *     .. Parameters ..
343        REAL ONE , ZERO
344        PARAMETER( ONE = 1.0E + 0 )
345        PARAMETER( ZERO = 0.0E + 0 )
346        COMPLEX CONE , CZERO
347        PARAMETER( CONE =( 1.0E + 0 , 0.0E + 0 ) )
348        PARAMETER( CZERO =( 0.0E + 0 , 0.0E + 0 ) )
349        INTEGER INT_ONE
350        PARAMETER( INT_ONE = 1 )
351        INTEGER DESCMULT , BIGNUM
352        PARAMETER(DESCMULT = 100 , BIGNUM = DESCMULT * DESCMULT)
353        INTEGER BLOCK_CYCLIC_2D , CSRC_ , CTXT_ , DLEN_ , DTYPE_ ,
354       $LLD_ , MB_ , M_ , NB_ , N_ , RSRC_
355        PARAMETER( BLOCK_CYCLIC_2D = 1 , DLEN_ = 9 , DTYPE_ = 1 ,
356       $CTXT_ = 2 , M_ = 3 , N_ = 4 , MB_ = 5 , NB_ = 6 ,
357       $RSRC_ = 7 , CSRC_ = 8 , LLD_ = 9 )
358        IF( INFO.EQ.0 ) THEN
359  
360  *         Use diagonal block(s) to modify this offdiagonal block
361  
362            CALL CTRSM( 'R' , 'L' , 'C' , 'N' , BW , BW , CONE ,
363       $    AF( ODD_SIZE*BW + MBW2 + 1 ) , BW ,
364       $    AF( ODD_SIZE*BW + 2*MBW2 + 1 ) , BW )
365  
366        ENDIF
367  *     End of "if( info.eq.0 ) then"
368  
369  *     Use offdiag block(s) to calculate modification to diag block
370  *     of processor to the left
371  
372        CALL CHERK( 'L' , 'N' , BW , BW , - ONE ,
373       $AF((ODD_SIZE + 2*BW)*BW + 1 ) , BW , ZERO ,
374       $WORK( 1 ) , BW )
375  
376  *     Send contribution to diagonal block's owning processor.
377  
378        CALL CGESD2D( ICTXT , BW , BW , WORK( 1 ) , BW ,
379       $0 , MYCOL - LEVEL_DIST )
380  
381  *     *******************************************************
382  
383        IF( MYCOL / LEVEL_DIST .LE.(NPCOL - 1) / LEVEL_DIST - 2 ) THEN
384  
385  *         Decide which processor offdiagonal block(s) goes to
386  
387            IF(( MOD( MYCOL / ( 2*LEVEL_DIST ) , 2 )) .EQ.0 ) THEN
388                COMM_PROC = MYCOL + LEVEL_DIST
389            ELSE
390                COMM_PROC = MYCOL - LEVEL_DIST
391            ENDIF
392  
393  *         Use offdiagonal blocks to calculate offdiag
394  *         block to send to neighboring processor. Depending
395  *         on circumstances , may need to transpose the matrix.
396  
397            CALL CGEMM( 'N' , 'N' , BW , BW , BW , - CONE ,
398       $    AF( ODD_SIZE*BW + 2*MBW2 + 1 ) , BW ,
399       $    AF( ODD_SIZE*BW + 1 ) , BW , CZERO , WORK( 1 ) ,
400       $    BW )
401  
402  *         Send contribution to offdiagonal block's owning processor.
403  
404            CALL CGESD2D( ICTXT , BW , BW , WORK( 1 ) , BW ,
405       $    0 , COMM_PROC )
406  
407        ENDIF
408  
409        ENDIF
410  *     End of "if( mycol/level_dist.le.(npcol-1)/level_dist -1 )..."
411  
412     14 CONTINUE
413  
414        ELSE
415  
416  *         CASE UPLO = 'U'
417  
418  *         *******************************************************************
419  *         PHASE 1 : Local computation phase.
420  *         *******************************************************************
421  
422  *         Sizes of the extra triangles communicated bewtween processors
423  
424            IF( MYCOL .GT. 0 ) THEN
425                PREV_TRI_SIZE_M = MIN( BW ,
426       $        NUMROC( N , PART_SIZE , MYCOL , 0 , NPCOL ) )
427                PREV_TRI_SIZE_N = MIN( BW ,
428       $        NUMROC( N , PART_SIZE , MYCOL - 1 , 0 , NPCOL ) )
429            ENDIF
430  
431            IF( MYCOL .LT. NPCOL - 1 ) THEN
432                NEXT_TRI_SIZE_M = MIN( BW ,
433       $        NUMROC( N , PART_SIZE , MYCOL + 1 , 0 , NPCOL ) )
434                NEXT_TRI_SIZE_N = MIN( BW ,
435       $        NUMROC( N , PART_SIZE , MYCOL , 0 , NPCOL ) )
436            ENDIF
437  
438  *         Factor main partition A_i^C = U_i {U_i}^C in each processor
439  
440            CALL CPBTRF( UPLO , ODD_SIZE , BW , A( OFST + 1) ,
441       $    LLDA , INFO )
442  
443            IF( INFO.NE.0 ) THEN
444                INFO = MYCOL + 1
445                GOTO 1600
446            ENDIF
447  
448            IF( MYCOL .LT. NP - 1 ) THEN
449  *             Apply factorization to odd - even connection block B_i
450  
451  *             Move the connection block in preparation.
452  
453                CALL CLACPY( 'L' , BW , BW , A(( OFST + 1 + ODD_SIZE*LLDA ) ) ,
454       $        LLDA - 1 , AF( ODD_SIZE*BW + 2*MBW2 + 1 + BW - BW ) , BW )
455  
456  *             Perform the triangular solve {L_i}{{B'}_i}^C = {B_i}^C
457  
458                CALL CTRTRS( 'U' , 'C' , 'N' , BW , BW ,
459       $        A( OFST + BW + 1 + (ODD_SIZE - BW)*LLDA ) , LLDA - 1 ,
460       $        AF( ODD_SIZE*BW + 2*MBW2 + 1 ) , BW , INFO )
461  
462  *             Move the resulting block back to its location in main storage.
463  
464                CALL CLACPY( 'L' , BW , BW , AF( ODD_SIZE*BW + 2*MBW2 + 1 + BW - BW ) ,
465       $        BW , A(( OFST + 1 + ODD_SIZE*LLDA )) , LLDA - 1 )
466  
467  *             Compute contribution to diagonal block(s) of reduced system.
468  *             {C'}_i^C = {C_i}^C - {{B'}_i}^C{{B'}_i}
469  
470  *             The following method uses more flops than necessary but
471  *             does not necessitate the writing of a new BLAS routine.
472  
473                CALL CHERK( UPLO , 'C' , BW , BW , - ONE ,
474       $        AF( ODD_SIZE*BW + 2*MBW2 + 1 ) , BW , ONE ,
475       $        A( OFST + BW + 1 + ODD_SIZE*LLDA ) , LLDA - 1 )
476  
477            ENDIF
478  *         End of "if( MYCOL .lt. NP-1 )..." loop
479  
480   1600 CONTINUE
481  *     If the processor could not locally factor , it jumps here.
482  
483        IF( MYCOL .NE. 0 ) THEN
484  *         Discard temporary matrix stored beginning in
485  *         AF((odd_size + 2*bw)*bw + 1 ) and use for
486  *         off_diagonal block of reduced system.
487  
488  *         Calculate the "spike" fillin , ${L_i} {{G}_i}^C = {D_i}$ .
489  
490  *         Copy D block into AF storage for solve.
491  
492            CALL CLATCPY( 'L' , PREV_TRI_SIZE_N , PREV_TRI_SIZE_M ,
493       $    A( OFST + 1 ) , LLDA - 1 , AF( 1 ) , ODD_SIZE )
494  
495            IF( INFO.EQ.0 ) THEN
496  
497                CALL CTBTRS( 'U' , 'C' , 'N' , ODD_SIZE , BW , BW ,
498       $        A( OFST + 1 ) , LLDA ,
499       $        AF( 1 ) , ODD_SIZE , INFO )
500  
501  *             Calculate the update block for previous proc , E_i = G_i{G_i}^C
502  
503                CALL CHERK( 'L' , 'C' , BW , ODD_SIZE ,
504       $        - ONE , AF( 1 ) , ODD_SIZE , ZERO ,
505       $        AF( 1 + (ODD_SIZE + 2*BW)*BW) , BW )
506  
507  *             Initiate send of E_i to previous processor to overlap
508  *             with next computation.
509  
510                CALL CGESD2D( ICTXT , BW , BW , AF( ODD_SIZE*BW + 2*MBW2 + 1 ) , BW ,
511       $        0 , MYCOL - 1 )
512  
513                IF( MYCOL .LT. NP - 1 ) THEN
514  
515  *                 Calculate off - diagonal block(s) of reduced system.
516  *                 Note : for ease of use in solution of reduced system , store
517  *                 L's off - diagonal block in conjugate transpose form.
518  *                 {F_i}^C = {H_i}{{B'}_i}^C
519  
520  *                 Copy matrix H_i(the last bw cols of G_i) to AF storage
521  *                 as per requirements of BLAS routine CTRMM.
522  *                 Since we have G_i^C stored , conjugate transpose
523  *                 H_i^C to H_i.
524  
525                    CALL CLATCPY( 'N' , BW , BW ,
526       $            AF( ODD_SIZE - BW + 1 ) , ODD_SIZE ,
527       $            AF((ODD_SIZE)*BW + 1) , BW )
528  
529                    CALL CTRMM( 'R' , 'L' , 'N' , 'N' , BW , BW , - CONE ,
530       $            A(( OFST + 1 + ODD_SIZE*LLDA ) ) , LLDA - 1 ,
531       $            AF((ODD_SIZE)*BW + 1 ) , BW )
532  
533                ENDIF
534  
535            ENDIF
536  *         End of "if( MYCOL .ne. 0 )..."
537  
538        ENDIF
539  *     End of "if(info.eq.0) then"
540  
541  *     Check to make sure no processors have found errors
542  
543        CALL IGAMX2D( ICTXT , 'A' , ' ' , 1 , 1 , INFO , 1 , INFO , INFO ,
544       $- 1 , 0 , 0 )
545  
546        IF( MYCOL.EQ.0 ) THEN
547            CALL IGEBS2D( ICTXT , 'A' , ' ' , 1 , 1 , INFO , 1 )
548        ELSE
549            CALL IGEBR2D( ICTXT , 'A' , ' ' , 1 , 1 , INFO , 1 , 0 , 0 )
550        ENDIF
551  
552        IF( INFO.NE.0 ) THEN
553            GOTO 1000
554        ENDIF
555  *     No errors found , continue
556  
557  *     *******************************************************************
558  *     PHASE 2 : Formation and factorization of Reduced System.
559  *     *******************************************************************
560  
561  *     Gather up local sections of reduced system
562  
563  *     The last processor does not participate in the factorization of
564  *     the reduced system , having sent its E_i already.
565        IF( MYCOL .EQ. NPCOL - 1 ) THEN
566            GOTO 24
567        ENDIF
568  
569  *     Initiate send of off - diag block(s) to overlap with next part.
570  *     Off - diagonal block needed on neighboring processor to start
571  *     algorithm.
572  
573        IF((MOD( MYCOL + 1 , 2 ) .EQ. 0) .AND.( MYCOL .GT. 0 ) ) THEN
574  
575            CALL CGESD2D( ICTXT , BW , BW ,
576       $    AF( ODD_SIZE*BW + 1 ) ,
577       $    BW , 0 , MYCOL - 1 )
578  
579        ENDIF
580  
581  *     Transpose last diagonal block into AF storage for subsequent
582  *     operations.
583  
584        CALL CLATCPY( 'U' , BW , BW ,
585       $A( OFST + ODD_SIZE*LLDA + 1 + BW ) ,
586       $LLDA - 1 , AF( ODD_SIZE*BW + MBW2 + 1 ) ,
587       $BW )
588  
589  *     Receive cont. to diagonal block that is stored on this proc.
590  
591        IF( MYCOL.LT. NPCOL - 1 ) THEN
592  
593            CALL CGERV2D( ICTXT , BW , BW ,
594       $    AF( ODD_SIZE*BW + 2*MBW2 + 1 ) ,
595       $    BW , 0 ,
596       $    MYCOL + 1 )
597  
598  *         Add contribution to diagonal block
599  
600            CALL CAXPY( MBW2 , CONE ,
601       $    AF( ODD_SIZE*BW + 2*MBW2 + 1 ) ,
602       $    1 , AF( ODD_SIZE*BW + MBW2 + 1 ) , 1 )
603  
604        ENDIF
605  
606  *     *************************************
607  *     Modification Loop
608  
609  *     The distance for sending and receiving for each level starts
610  *     at 1 for the first level.
611        LEVEL_DIST = 1
612  
613  *     Do until this proc is needed to modify other procs' equations
614  
615     22 CONTINUE
616        IF( MOD((MYCOL + 1) / LEVEL_DIST , 2) .NE. 0 ) GOTO 21
617  
618  *     Receive and add contribution to diagonal block from the left
619  
620        IF( MYCOL - LEVEL_DIST .GE. 0 ) THEN
621            CALL CGERV2D( ICTXT , BW , BW , WORK( 1 ) ,
622       $    BW , 0 , MYCOL - LEVEL_DIST )
623  
624            CALL CAXPY( MBW2 , CONE , WORK( 1 ) , 1 ,
625       $    AF( ODD_SIZE*BW + MBW2 + 1 ) , 1 )
626  
627        ENDIF
628  
629  *     Receive and add contribution to diagonal block from the right
630  
631        IF( MYCOL + LEVEL_DIST .LT. NPCOL - 1 ) THEN
632            CALL CGERV2D( ICTXT , BW , BW , WORK( 1 ) ,
633       $    BW , 0 , MYCOL + LEVEL_DIST )
634  
635            CALL CAXPY( MBW2 , CONE , WORK( 1 ) , 1 ,
636       $    AF( ODD_SIZE*BW + MBW2 + 1 ) , 1 )
637  
638        ENDIF
639  
640        LEVEL_DIST = LEVEL_DIST*2
641  
642        GOTO 22
643     21 CONTINUE
644  *     [End of GOTO Loop]
645  
646  *     *********************************
647  *     Calculate and use this proc's blocks to modify other procs'...
648  
649  *     Factor diagonal block
650  
651        CALL CPOTRF( 'L' , BW , AF( ODD_SIZE*BW + MBW2 + 1 ) ,
652       $BW , INFO )
653  
654        IF( INFO.NE.0 ) THEN
655            INFO = NPCOL + MYCOL
656        ENDIF
657  
658  *     ****************************************************************
659  *     Receive offdiagonal block from processor to right.
660  *     If this is the first group of processors , the receive comes
661  *     from a different processor than otherwise.
662  
663        IF( LEVEL_DIST .EQ. 1 )THEN
664            COMM_PROC = MYCOL + 1
665  
666  *         Move block into place that it will be expected to be for
667  *         calcs.
668  
669            CALL CLACPY( 'N' , BW , BW , AF( ODD_SIZE*BW + 1 ) , BW ,
670       $    AF( ODD_SIZE*BW + 2*MBW2 + 1 ) , BW )
671  
672        ELSE
673            COMM_PROC = MYCOL + LEVEL_DIST / 2
674        ENDIF
675  
676        IF( MYCOL / LEVEL_DIST .LE.(NPCOL - 1) / LEVEL_DIST - 2 )THEN
677  
678            CALL CGERV2D( ICTXT , BW , BW ,
679       $    AF( ODD_SIZE*BW + 1 ) ,
680       $    BW , 0 , COMM_PROC )
681  
682            IF( INFO .EQ. 0 ) THEN
683  
684  *             Modify upper off_diagonal block with diagonal block
685  
686                CALL CTRSM( 'L' , 'L' , 'N' , 'N' , BW , BW , CONE ,
687       $        AF( ODD_SIZE*BW + MBW2 + 1 ) , BW ,
688       $        AF( ODD_SIZE*BW + 1 ) , BW )
689  
690            ENDIF
691  *         End of "if( info.eq.0 ) then"
692  
693  *         Calculate contribution from this block to next diagonal block
694  
695            CALL CHERK( 'L' , 'C' , BW , BW , - ONE ,
696       $    AF((ODD_SIZE)*BW + 1 ) , BW , ZERO ,
697       $    WORK( 1 ) , BW )
698  
699  *         Send contribution to diagonal block's owning processor.
700  
701            CALL CGESD2D( ICTXT , BW , BW , WORK( 1 ) , BW ,
702       $    0 , MYCOL + LEVEL_DIST )
703  
704        ENDIF
705  *     End of "if( mycol/level_dist .le.(npcol-1)/level_dist-2 )..."
706  
707  *     ****************************************************************
708  *     Receive off_diagonal block from left and use to finish with this
709  *     processor.
710  
711        IF((MYCOL / LEVEL_DIST .GT. 0 ).AND.
712       $( MYCOL / LEVEL_DIST .LE.(NPCOL - 1) / LEVEL_DIST - 1 ) ) THEN
713  
714        IF( LEVEL_DIST .GT. 1)THEN
715  
716  *         Receive offdiagonal block(s) from proc level_dist / 2 to the
717  *         left
718  
719            CALL CGERV2D( ICTXT , BW , BW ,
720       $    AF( ODD_SIZE*BW + 2*MBW2 + 1 ) ,
721       $    BW , 0 , MYCOL - LEVEL_DIST / 2 )
722  
723        ENDIF
724  
725        IF( INFO.EQ.0 ) THEN
726  
727  *         Use diagonal block(s) to modify this offdiagonal block
728  
729            CALL CTRSM( 'R' , 'L' , 'C' , 'N' , BW , BW , CONE ,
730       $    AF( ODD_SIZE*BW + MBW2 + 1 ) , BW ,
731       $    AF( ODD_SIZE*BW + 2*MBW2 + 1 ) , BW )
732  
733        ENDIF
734  *     End of "if( info.eq.0 ) then"
735  
736  *     Use offdiag block(s) to calculate modification to diag block
737  *     of processor to the left
738  
739        CALL CHERK( 'L' , 'N' , BW , BW , - ONE ,
740       $AF((ODD_SIZE + 2*BW)*BW + 1 ) , BW , ZERO ,
741       $WORK( 1 ) , BW )
742  
743  *     Send contribution to diagonal block's owning processor.
744  
745        CALL CGESD2D( ICTXT , BW , BW , WORK( 1 ) , BW ,
746       $0 , MYCOL - LEVEL_DIST )
747  
748  *     *******************************************************
749  
750        IF( MYCOL / LEVEL_DIST .LE.(NPCOL - 1) / LEVEL_DIST - 2 ) THEN
751  
752  *         Decide which processor offdiagonal block(s) goes to
753  
754            IF(( MOD( MYCOL / ( 2*LEVEL_DIST ) , 2 )) .EQ.0 ) THEN
755                COMM_PROC = MYCOL + LEVEL_DIST
756            ELSE
757                COMM_PROC = MYCOL - LEVEL_DIST
758            ENDIF
759  
760  *         Use offdiagonal blocks to calculate offdiag
761  *         block to send to neighboring processor. Depending
762  *         on circumstances , may need to transpose the matrix.
763  
764            CALL CGEMM( 'N' , 'N' , BW , BW , BW , - CONE ,
765       $    AF( ODD_SIZE*BW + 2*MBW2 + 1 ) , BW ,
766       $    AF( ODD_SIZE*BW + 1 ) , BW , CZERO , WORK( 1 ) ,
767       $    BW )
768  
769  *         Send contribution to offdiagonal block's owning processor.
770  
771            CALL CGESD2D( ICTXT , BW , BW , WORK( 1 ) , BW ,
772       $    0 , COMM_PROC )
773  
774        ENDIF
775  
776        ENDIF
777  *     End of "if( mycol/level_dist.le.(npcol-1)/level_dist -1 )..."
778  
779     24 CONTINUE
780  
781        ENDIF
782  
783   1000 CONTINUE
784  
785  *     Free BLACS space used to hold standard - form grid.
786  
787        IF( ICTXT_SAVE .NE. ICTXT_NEW ) THEN
788            CALL BLACS_GRIDEXIT( ICTXT_NEW )
789        ENDIF
790  
791   1234 CONTINUE
792  
793  *     Restore saved input parameters
794  
795        ICTXT = ICTXT_SAVE
796        NP = NP_SAVE
797  
798  *     Output minimum worksize
799  
800        WORK( 1 ) = WORK_SIZE_MIN
801  
802  *     Make INFO consistent across processors
803  
804        CALL IGAMX2D( ICTXT , 'A' , ' ' , 1 , 1 , INFO , 1 , INFO , INFO ,
805       $- 1 , 0 , 0 )
806  
807        IF( MYCOL.EQ.0 ) THEN
808            CALL IGEBS2D( ICTXT , 'A' , ' ' , 1 , 1 , INFO , 1 )
809        ELSE
810            CALL IGEBR2D( ICTXT , 'A' , ' ' , 1 , 1 , INFO , 1 , 0 , 0 )
811        ENDIF
812  
813        RETURN
814  
815  *     End of PCPBTRF
816  
817        END