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