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