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