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