001        SUBROUTINE PSPBTRF( 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        REAL A( * ) , AF( * ) , WORK( * )
016  *     ..
017  
018  *     Purpose
019  *     === ====
020  
021  *     PSPBTRF 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 PSPBTRS 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) REAL 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) REAL array , dimension LAF.
085  *     Auxiliary Fillin Space.
086  *     Fillin is created during the factorization routine
087  *     PSPBTRF and this is stored in AF. If a linear system
088  *     is to be solved using PSPBTRS 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  *     REAL 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        REAL ONE
342        PARAMETER( ONE = 1.0E + 0 )
343        REAL ZERO
344        PARAMETER( ZERO = 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        END IF
355  
356        IF( INFO.EQ.0 ) THEN
357  
358  *         Use diagonal block(s) to modify this offdiagonal block
359  
359               
360            CALL STRSM( '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 SSYRK( '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 SGESD2D( 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  
384               
385            IF(( MOD( MYCOL / ( 2*LEVEL_DIST ) , 2 ) ).EQ.0 ) THEN
385                   
386                COMM_PROC = MYCOL + LEVEL_DIST
387            ELSE
387                   
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 SGEMM( '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 SGESD2D( 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  
421               
422            IF( MYCOL.GT.0 ) THEN
422                   
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
429                   
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 SPBTRF( UPLO , ODD_SIZE , BW , A( OFST + 1 ) , LLDA , INFO )
439  
440            IF( INFO.NE.0 ) THEN
440                   
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  
449                   
450                CALL SLACPY( '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 STRTRS( '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 SLACPY( '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 SSYRK( 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  
488               
489            CALL SLATCPY( '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  
493                   
494                CALL STBTRS( '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 SSYRK( '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 SGESD2D( 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 STRMM.
518  *                 Since we have G_i^T stored , transpose
519  *                 H_i^T to H_i.
520  
520                       
521                    CALL SLATCPY( 'N' , BW , BW , AF( ODD_SIZE - BW + 1 ) ,
522       $            ODD_SIZE , AF(( ODD_SIZE )*BW + 1 ) , BW )
523  
524                    CALL STRMM( '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
541               
542            CALL IGEBS2D( ICTXT , 'A' , ' ' , 1 , 1 , INFO , 1 )
543        ELSE
543               
544            CALL IGEBR2D( ICTXT , 'A' , ' ' , 1 , 1 , INFO , 1 , 0 , 0 )
545        END IF
546  
547        IF( INFO.NE.0 ) THEN
547               
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
560               
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  
569               
570            CALL SGESD2D( 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 SLATCPY( '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  
584               
585            CALL SGERV2D( 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 SAXPY( 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 )
605               
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
610                   
611                CALL SGERV2D( ICTXT , BW , BW , WORK( 1 ) , BW , 0 ,
612       $        MYCOL - LEVEL_DIST )
613  
614                CALL SAXPY( 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
621                   
622                CALL SGERV2D( ICTXT , BW , BW , WORK( 1 ) , BW , 0 ,
623       $        MYCOL + LEVEL_DIST )
624  
625                CALL SAXPY( 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 SPOTRF( 'L' , BW , AF( ODD_SIZE*BW + MBW2 + 1 ) , BW , INFO )
642  
643        IF( INFO.NE.0 ) THEN
643               
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
652               
653            COMM_PROC = MYCOL + 1
654  
655  *         Move block into place that it will be expected to be for
656  *         calcs.
657  
658            CALL SLACPY( 'N' , BW , BW , AF( ODD_SIZE*BW + 1 ) , BW ,
659       $    AF( ODD_SIZE*BW + 2*MBW2 + 1 ) , BW )
660  
661        ELSE
661               
662            COMM_PROC = MYCOL + LEVEL_DIST / 2
663        END IF
664  
665        IF( MYCOL / LEVEL_DIST.LE.( NPCOL - 1 ) / LEVEL_DIST - 2 ) THEN
666  
666               
667            CALL SGERV2D( 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  
673                   
674                CALL STRSM( '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 SSYRK( '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 SGESD2D( 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  
705               
706            CALL SGERV2D( 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  
714               
715            CALL STRSM( '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 SSYRK( '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 SGESD2D( 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  
739               
740            IF(( MOD( MYCOL / ( 2*LEVEL_DIST ) , 2 ) ).EQ.0 ) THEN
740                   
741                COMM_PROC = MYCOL + LEVEL_DIST
742            ELSE
742                   
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 SGEMM( '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 SGESD2D( 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
773               
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
793               
794            CALL IGEBS2D( ICTXT , 'A' , ' ' , 1 , 1 , INFO , 1 )
795        ELSE
795               
796            CALL IGEBR2D( ICTXT , 'A' , ' ' , 1 , 1 , INFO , 1 , 0 , 0 )
797        END IF
798  
799        RETURN
800  
801  *     End of PSPBTRF
802  
803        END387203