Routine: PSPTTRF()  File: SRC\pspttrf.f

 
 
# lines: 1017
  # code: 1017
  # comment: 0
  # blank:0
# Variables:28
# Callers:1
# Callings:0
# Words:54
# Keywords:32
 

 

 
..
     .. Local Scalars ..
     ..
     .. Local Arrays ..
     ..
     .. External Subroutines ..
     ..
     .. External Functions ..
     ..
     .. Intrinsic Functions ..
     ..
     .. Executable Statements ..
     Test the input parameters
     Convert descriptor into standard form for easy access to
        parameters, check that grid is of right shape.
        Temporarily set the descriptor type to 1xP type
     Get values out of descriptor for use in code.
     Get grid parameters
     Argument checking that is specific to Divide & Conquer routine
     Check auxiliary storage size
        put minimum value of laf into AF( 1 )
     Check worksize
     Pack params and positions into arrays for global consistency check
     Want to find errors with MIN( ), so if no error, set it to a big
     number. If there already is an error, multiply by the the
     descriptor multiplier.
     Check consistency across processors
     Prepare output: set info = 0 if no error, and divide by DESCMULT
     if error is not in a descriptor entry.
     Quick return if possible
     Adjust addressing into matrix space to properly get into
        the beginning part of the relevant data
     Form a new BLACS grid (the "standard form" grid) with only procs
        holding part of the matrix, of size 1xNP where NP is adjusted,
        starting at csrc=0, with JA modified to reflect dropped procs.
     First processor to hold part of the matrix:
     Calculate new JA one while dropping off unused processors.
     Save and compute new value of NP
     Call utility routine that forms "standard-form" grid
     Use new context from standard grid as context.
     Get information about new grid.
     Drop out processors that do not have part of the matrix.
     ********************************
     Values reused throughout routine
     User-input value of partition size
     Number of columns in each processor
     Offset in columns to beginning of main partition in each proc
     Size of main (or odd) partition in each processor
       Zero out space for fillin
     Begin main code
*******************************************************************
       PHASE 1: Local computation phase.
*******************************************************************
         Transfer last triangle D_i of local matrix to next processor
         which needs it to calculate fillin due to factorization of
         its main (odd) block A_i.
         Overlap the send with the factorization of A_i.
       Factor main partition A_i = L_i {L_i}^T in each processor
       Or A_i = {U_i}^T {U_i} if E is the upper superdiagonal
         Apply factorization to odd-even connection block B_i
         Perform the triangular system solve {L_i}{{B'}_i}^T = {B_i}^T
           by dividing B_i by diagonal element
         Compute contribution to diagonal block(s) of reduced system.
          {C'}_i = {C_i}-{{B'}_i}{{B'}_i}^T
       End of "if ( MYCOL .lt. NP-1 )..." loop
       If the processor could not locally factor, it jumps here.
         Receive previously transmitted matrix section, which forms
         the right-hand-side for the triangular solve that calculates
         the "spike" fillin.
         Calculate the "spike" fillin, ${L_i} {{G}_i}^T = {D_i}$ .
         Divide by D
         Calculate the update block for previous proc, E_i = G_i{G_i}^T
         Since there is no element-by-element vector multiplication in
           the BLAS, this loop must be hardwired in without a BLAS call
         Initiate send of E_i to previous processor to overlap
           with next computation.
           Calculate off-diagonal block(s) of reduced system.
           Note: for ease of use in solution of reduced system, store
           L's off-diagonal block in transpose form.
           {F_i}^T =  {H_i}{{B'}_i}^T
       End of "if ( MYCOL .ne. 0 )..."
       End of "if (info.eq.0) then"
       Check to make sure no processors have found errors
       No errors found, continue
*******************************************************************
       PHASE 2: Formation and factorization of Reduced System.
*******************************************************************
       Gather up local sections of reduced system
     The last processor does not participate in the factorization of
       the reduced system, having sent its E_i already.
       Initiate send of off-diag block(s) to overlap with next part.
       Off-diagonal block needed on neighboring processor to start
       algorithm.
       Copy last diagonal block into AF storage for subsequent
         operations.
       Receive cont. to diagonal block that is stored on this proc.
          Add contribution to diagonal block
       *************************************
       Modification Loop
       The distance for sending and receiving for each level starts
         at 1 for the first level.
       Do until this proc is needed to modify other procs' equations
         Receive and add contribution to diagonal block from the left
         Receive and add contribution to diagonal block from the right
       [End of GOTO Loop]
       *********************************
       Calculate and use this proc's blocks to modify other procs'...
       ****************************************************************
       Receive offdiagonal block from processor to right.
         If this is the first group of processors, the receive comes
         from a different processor than otherwise.
           Move block into place that it will be expected to be for
             calcs.
         Modify upper off_diagonal block with diagonal block
         End of "if ( info.eq.0 ) then"
         Calculate contribution from this block to next diagonal block
         Send contribution to diagonal block's owning processor.
       End of "if( mycol/level_dist .le. (npcol-1)/level_dist-2 )..."
       ****************************************************************
       Receive off_diagonal block from left and use to finish with this
         processor.
           Receive offdiagonal block(s) from proc level_dist/2 to the
           left
         Use diagonal block(s) to modify this offdiagonal block
         End of "if( info.eq.0 ) then"
         Use offdiag block(s) to calculate modification to diag block
           of processor to the left
         Send contribution to diagonal block's owning processor.
         *******************************************************
           Decide which processor offdiagonal block(s) goes to
           Use offdiagonal blocks to calculate offdiag
             block to send to neighboring processor. Depending
             on circumstances, may need to transpose the matrix.
           Send contribution to offdiagonal block's owning processor.
 
Display dynamic version Find AutoScroll Reload FontSize: - + Hide Comments Hide Blanks Frame FullScreen MailPrint

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