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