001        SUBROUTINE PDPTTRSV( UPLO , N , NRHS , D , E , JA , DESCA , B , IB , DESCB ,
002       $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 UPLO
011        INTEGER IB , INFO , JA , LAF , LWORK , N , NRHS
012  *     ..
013  *     .. Array Arguments ..
014        INTEGER DESCA( * ) , DESCB( * )
015        DOUBLE PRECISION AF( * ) , B( * ) , D( * ) , E( * ) , WORK( * )
016  *     ..
017  
018  *     Purpose
019  *     === ====
020  
021  *     PDPTTRSV solves a tridiagonal 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)^T * X = B(IB : IB + N - 1 , 1 : NRHS)
026  
027  *     where A(1 : N , JA : JA + N - 1) is a tridiagonal
028  *     triangular matrix factor produced by the
029  *     Cholesky factorization code PDPTTRF
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)^T
034  *     is dictated by the user by the parameter TRANS.
035  
036  *     Routine PDPTTRF 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  *     N(global input) INTEGER
048  *     The number of rows and columns to be operated on , i.e. the
049  *     order of the distributed submatrix A(1 : N , JA : JA + N - 1). N >= 0.
050  
051  *     NRHS(global input) INTEGER
052  *     The number of right hand sides , i.e. , the number of columns
053  *     of the distributed submatrix B(IB : IB + N - 1 , 1 : NRHS).
054  *     NRHS >= 0.
055  
056  *     D(local input / local output) DOUBLE PRECISION pointer to local
057  *     part of global vector storing the main diagonal of the
058  *     matrix.
059  *     On exit , this array contains information containing the
060  *     factors of the matrix.
061  *     Must be of size >= DESCA( NB_ ).
062  
063  *     E(local input / local output) DOUBLE PRECISION pointer to local
064  *     part of global vector storing the upper diagonal of the
065  *     matrix. Globally , DU(n) is not referenced , and DU must be
066  *     aligned with D.
067  *     On exit , this array contains information containing the
068  *     factors of the matrix.
069  *     Must be of size >= DESCA( NB_ ).
070  
071  *     JA(global input) INTEGER
072  *     The index in the global array A that points to the start of
073  *     the matrix to be operated on(which may be either all of A
074  *     or a submatrix of A).
075  
076  *     DESCA(global and local input) INTEGER array of dimension DLEN.
077  *     if 1D type(DTYPE_A = 501 or 502) , DLEN >= 7 ;
078  *     if 2D type(DTYPE_A = 1) , DLEN >= 9.
079  *     The array descriptor for the distributed matrix A.
080  *     Contains information of mapping of A to memory. Please
081  *     see NOTES below for full description and options.
082  
083  *     B(local input / local output) DOUBLE PRECISION pointer into
084  *     local memory to an array of local lead dimension lld_b >= NB.
085  *     On entry , this array contains the
086  *     the local pieces of the right hand sides
087  *     B(IB : IB + N - 1 , 1 : NRHS).
088  *     On exit , this contains the local piece of the solutions
089  *     distributed matrix X.
090  
091  *     IB(global input) INTEGER
092  *     The row index in the global array B that points to the first
093  *     row of the matrix to be operated on(which may be either
094  *     all of B or a submatrix of B).
095  
096  *     DESCB(global and local input) INTEGER array of dimension DLEN.
097  *     if 1D type(DTYPE_B = 502) , DLEN >= 7 ;
098  *     if 2D type(DTYPE_B = 1) , DLEN >= 9.
099  *     The array descriptor for the distributed matrix B.
100  *     Contains information of mapping of B to memory. Please
101  *     see NOTES below for full description and options.
102  
103  *     AF(local output) DOUBLE PRECISION array , dimension LAF.
104  *     Auxiliary Fillin Space.
105  *     Fillin is created during the factorization routine
106  *     PDPTTRF and this is stored in AF. If a linear system
107  *     is to be solved using PDPTTRS after the factorization
108  *     routine , AF *must not be altered* after the factorization.
109  
110  *     LAF(local input) INTEGER
111  *     Size of user - input Auxiliary Fillin space AF. Must be >=
112  *     (NB + 2)
113  *     If LAF is not large enough , an error code will be returned
114  *     and the minimum acceptable size will be returned in AF( 1 )
115  
116  *     WORK(local workspace / local output)
117  *     DOUBLE PRECISION temporary workspace. This space may
118  *     be overwritten in between calls to routines. WORK must be
119  *     the size given in LWORK.
120  *     On exit , WORK( 1 ) contains the minimal LWORK.
121  
122  *     LWORK(local input or global input) INTEGER
123  *     Size of user - input workspace WORK.
124  *     If LWORK is too small , the minimal acceptable size will be
125  *     returned in WORK(1) and an error code is returned. LWORK >=
126  *     (10 + 2*min(100 , NRHS))*NPCOL + 4*NRHS
127  
128  *     INFO(local output) INTEGER
129  *     = 0 : successful exit
130  *     < 0 : If the i - th argument is an array and the j - entry had
131  *     an illegal value , then INFO = - (i*100 + j) , if the i - th
132  *     argument is a scalar and had an illegal value , then
133  *     INFO = - i.
134  
135  *     === ==================================================================
136  
137  *     Restrictions
138  *     === =========
139  
140  *     The following are restrictions on the input parameters. Some of these
141  *     are temporary and will be removed in future releases , while others
142  *     may reflect fundamental technical limitations.
143  
144  *     Non - cyclic restriction : VERY IMPORTANT !
145  *     P*NB >= mod(JA - 1 , NB) + N.
146  *     The mapping for matrices must be blocked , reflecting the nature
147  *     of the divide and conquer algorithm as a task - parallel algorithm.
148  *     This formula in words is : no processor may have more than one
149  *     chunk of the matrix.
150  
151  *     Blocksize cannot be too small :
152  *     If the matrix spans more than one processor , the following
153  *     restriction on NB , the size of each block on each processor ,
154  *     must hold :
155  *     NB >= 2
156  *     The bulk of parallel computation is done on the matrix of size
157  *     O(NB) on each processor. If this is too small , divide and conquer
158  *     is a poor choice of algorithm.
159  
160  *     Submatrix reference :
161  *     JA = IB
162  *     Alignment restriction that prevents unnecessary communication.
163  
164  *     === ==================================================================
165  
166  *     Notes
167  *     === ==
168  
169  *     If the factorization routine and the solve routine are to be called
170  *     separately(to solve various sets of righthand sides using the same
171  *     coefficient matrix) , the auxiliary space AF *must not be altered*
172  *     between calls to the factorization routine and the solve routine.
173  
174  *     The best algorithm for solving banded and tridiagonal linear systems
175  *     depends on a variety of parameters , especially the bandwidth.
176  *     Currently , only algorithms designed for the case N / P >> bw are
177  *     implemented. These go by many names , including Divide and Conquer ,
178  *     Partitioning , domain decomposition - type , etc.
179  *     For tridiagonal matrices , it is obvious : N / P >> bw( = 1) , and so D&C
180  *     algorithms are the appropriate choice.
181  
182  *     Algorithm description : Divide and Conquer
183  
184  *     The Divide and Conqer algorithm assumes the matrix is narrowly
185  *     banded compared with the number of equations. In this situation ,
186  *     it is best to distribute the input matrix A one - dimensionally ,
187  *     with columns atomic and rows divided amongst the processes.
188  *     The basic algorithm divides the tridiagonal matrix up into
189  *     P pieces with one stored on each processor ,
190  *     and then proceeds in 2 phases for the factorization or 3 for the
191  *     solution of a linear system.
192  *     1) Local Phase :
193  *     The individual pieces are factored independently and in
194  *     parallel. These factors are applied to the matrix creating
195  *     fillin , which is stored in a non - inspectable way in auxiliary
196  *     space AF. Mathematically , this is equivalent to reordering
197  *     the matrix A as P A P^T and then factoring the principal
198  *     leading submatrix of size equal to the sum of the sizes of
199  *     the matrices factored on each processor. The factors of
200  *     these submatrices overwrite the corresponding parts of A
201  *     in memory.
202  *     2) Reduced System Phase :
203  *     A small((P - 1)) system is formed representing
204  *     interaction of the larger blocks , and is stored(as are its
205  *     factors) in the space AF. A parallel Block Cyclic Reduction
206  *     algorithm is used. For a linear system , a parallel front solve
207  *     followed by an analagous backsolve , both using the structure
208  *     of the factored matrix , are performed.
209  *     3) Backsubsitution Phase :
210  *     For a linear system , a local backsubstitution is performed on
211  *     each processor in parallel.
212  
213  *     Descriptors
214  *     === ========
215  
216  *     Descriptors now have *types* and differ from ScaLAPACK 1.0.
217  
218  *     Note : tridiagonal codes can use either the old two dimensional
219  *     or new one - dimensional descriptors , though the processor grid in
220  *     both cases *must be one - dimensional*. We describe both types below.
221  
222  *     Each global data object is described by an associated description
223  *     vector. This vector stores the information required to establish
224  *     the mapping between an object element and its corresponding process
225  *     and memory location.
226  
227  *     Let A be a generic term for any 2D block cyclicly distributed array.
228  *     Such a global array has an associated description vector DESCA.
229  *     In the following comments , the character _ should be read as
230  *     "of the global array".
231  
232  *     NOTATION STORED IN EXPLANATION
233  *     --- ------------ -------------- --------------------------------------
234  *     DTYPE_A(global) DESCA( DTYPE_ )The descriptor type. In this case ,
235  *     DTYPE_A = 1.
236  *     CTXT_A(global) DESCA( CTXT_ ) The BLACS context handle , indicating
237  *     the BLACS process grid A is distribu -
238  *     ted over. The context itself is glo -
239  *     bal , but the handle(the integer
240  *     value) may vary.
241  *     M_A(global) DESCA( M_ ) The number of rows in the global
242  *     array A.
243  *     N_A(global) DESCA( N_ ) The number of columns in the global
244  *     array A.
245  *     MB_A(global) DESCA( MB_ ) The blocking factor used to distribute
246  *     the rows of the array.
247  *     NB_A(global) DESCA( NB_ ) The blocking factor used to distribute
248  *     the columns of the array.
249  *     RSRC_A(global) DESCA( RSRC_ ) The process row over which the first
250  *     row of the array A is distributed.
251  *     CSRC_A(global) DESCA( CSRC_ ) The process column over which the
252  *     first column of the array A is
253  *     distributed.
254  *     LLD_A(local) DESCA( LLD_ ) The leading dimension of the local
255  *     array. LLD_A >= MAX(1 , LOCr(M_A)).
256  
257  *     Let K be the number of rows or columns of a distributed matrix ,
258  *     and assume that its process grid has dimension p x q.
259  *     LOCr( K ) denotes the number of elements of K that a process
260  *     would receive if K were distributed over the p processes of its
261  *     process column.
262  *     Similarly , LOCc( K ) denotes the number of elements of K that a
263  *     process would receive if K were distributed over the q processes of
264  *     its process row.
265  *     The values of LOCr() and LOCc() may be determined via a call to the
266  *     ScaLAPACK tool function , NUMROC :
267  *     LOCr( M ) = NUMROC( M , MB_A , MYROW , RSRC_A , NPROW ) ,
268  *     LOCc( N ) = NUMROC( N , NB_A , MYCOL , CSRC_A , NPCOL ).
269  *     An upper bound for these quantities may be computed by :
270  *     LOCr( M ) <= ceil( ceil(M / MB_A) / NPROW )*MB_A
271  *     LOCc( N ) <= ceil( ceil(N / NB_A) / NPCOL )*NB_A
272  
273  *     One - dimensional descriptors :
274  
275  *     One - dimensional descriptors are a new addition to ScaLAPACK since
276  *     version 1.0. They simplify and shorten the descriptor for 1D
277  *     arrays.
278  
279  *     Since ScaLAPACK supports two - dimensional arrays as the fundamental
280  *     object , we allow 1D arrays to be distributed either over the
281  *     first dimension of the array(as if the grid were P - by - 1) or the
282  *     2nd dimension(as if the grid were 1 - by - P). This choice is
283  *     indicated by the descriptor type(501 or 502)
284  *     as described below.
285  *     However , for tridiagonal matrices , since the objects being
286  *     distributed are the individual vectors storing the diagonals , we
287  *     have adopted the convention that both the P - by - 1 descriptor and
288  *     the 1 - by - P descriptor are allowed and are equivalent for
289  *     tridiagonal matrices. Thus , for tridiagonal matrices ,
290  *     DTYPE_A = 501 or 502 can be used interchangeably
291  *     without any other change.
292  *     We require that the distributed vectors storing the diagonals of a
293  *     tridiagonal matrix be aligned with each other. Because of this , a
294  *     single descriptor , DESCA , serves to describe the distribution of
295  *     of all diagonals simultaneously.
296  
297  *     IMPORTANT NOTE : the actual BLACS grid represented by the
298  *     CTXT entry in the descriptor may be *either* P - by - 1 or 1 - by - P
299  *     irrespective of which one - dimensional descriptor type
300  *     (501 or 502) is input.
301  *     This routine will interpret the grid properly either way.
302  *     ScaLAPACK routines *do not support intercontext operations* so that
303  *     the grid passed to a single ScaLAPACK routine *must be the same*
304  *     for all array descriptors passed to that routine.
305  
306  *     NOTE : In all cases where 1D descriptors are used , 2D descriptors
307  *     may also be used , since a one - dimensional array is a special case
308  *     of a two - dimensional array with one dimension of size unity.
309  *     The two - dimensional array used in this case *must* be of the
310  *     proper orientation :
311  *     If the appropriate one - dimensional descriptor is DTYPEA = 501
312  *     (1 by P type) , then the two dimensional descriptor must
313  *     have a CTXT value that refers to a 1 by P BLACS grid ;
314  *     If the appropriate one - dimensional descriptor is DTYPEA = 502
315  *     (P by 1 type) , then the two dimensional descriptor must
316  *     have a CTXT value that refers to a P by 1 BLACS grid.
317  
318  *     Summary of allowed descriptors , types , and BLACS grids :
319  *     DTYPE 501 502 1 1
320  *     BLACS grid 1xP or Px1 1xP or Px1 1xP Px1
321  *     --- --------------------------------------------------
322  *     A               OK OK OK NO
323  *     B               NO OK NO OK
324  
325  *     Note that a consequence of this chart is that it is not possible
326  *     for *both* DTYPE_A and DTYPE_B to be 2D_type(1) , as these lead
327  *     to opposite requirements for the orientation of the BLACS grid ,
328  *     and as noted before , the *same* BLACS context must be used in
329  *     all descriptors in a single ScaLAPACK subroutine call.
330  
331  *     Let A be a generic term for any 1D block cyclicly distributed array.
332  *     Such a global array has an associated description vector DESCA.
333  *     In the following comments , the character _ should be read as
334  *     "of the global array".
335  
336  *     NOTATION STORED IN EXPLANATION
337  *     --- ------------ ---------- ------------------------------------------
338  *     DTYPE_A(global) DESCA( 1 ) The descriptor type. For 1D grids ,
339  *     TYPE_A = 501 : 1 - by - P grid.
340  *     TYPE_A = 502 : P - by - 1 grid.
341  *     CTXT_A(global) DESCA( 2 ) The BLACS context handle , indicating
342  *     the BLACS process grid A is distribu -
343  *     ted over. The context itself is glo -
344  *     bal , but the handle(the integer
345  *     value) may vary.
346  *     N_A(global) DESCA( 3 ) The size of the array dimension being
347  *     distributed.
348  *     NB_A(global) DESCA( 4 ) The blocking factor used to distribute
349  *     the distributed dimension of the array.
350  *     SRC_A(global) DESCA( 5 ) The process row or column over which the
351  *     first row or column of the array
352  *     is distributed.
353  *     Ignored DESCA( 6 ) Ignored for tridiagonal matrices.
354  *     Reserved DESCA( 7 ) Reserved for future use.
355  
356  *     === ==================================================================
357  
358  *     Code Developer : Andrew J. Cleary , University of Tennessee.
359  *     Current address : Lawrence Livermore National Labs.
360  
361  *     === ==================================================================
362  
363  *     .. Parameters ..
364        DOUBLE PRECISION ONE
365        PARAMETER( ONE = 1.0D + 0 )
366        DOUBLE PRECISION ZERO
367        PARAMETER( ZERO = 0.0D + 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        END30156