Routine: PSDBTRF()  File: SRC\psdbtrf.f

 
 
# lines: 1252
  # code: 1252
  # comment: 0
  # blank:0
# Variables:30
# Callers:1
# Callings:0
# Words:323
# Keywords:83
 

 

 
..
     .. 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.
     Get values out of descriptor for use in code.
     Get grid parameters
     Size of separator blocks is maximum of bandwidths
     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
     Offset in elements
     Size of main (or odd) partition in each processor
     Offset to workspace for Upper triangular factor
       Zero out space for fillin
       Zero out space for work
     Begin main code
*******************************************************************
       PHASE 1: Local computation phase.
*******************************************************************
       Sizes of the extra triangles communicated bewtween processors
         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 {U_i} in each processor
         Apply factorization to lower connection block BL_i
         transpose the connection block in preparation.
         Apply factorization to upper connection block BU_i
         Move the connection block in preparation.
         Perform the triangular system solve {L_i}{{BU'}_i} = {B_i}
         Perform the triangular solve {U_i}^T{BL'}_i^T = {BL_i}^T
         transpose resulting block to its location
           in main storage.
         Move the resulting block back to its location in main storage.
         Compute contribution to diagonal block(s) of reduced system.
          {C'}_i = {C_i}-{{BL'}_i}{{BU'}_i}
         The following method uses more flops than necessary but
           does not necessitate the writing of a new BLAS routine.
       End of "if ( MYCOL .lt. NP-1 )..." loop
       If the processor could not locally factor, it jumps here.
         Discard temporary matrix stored beginning in
           AF( (odd_size+2*bwl, bwu)*bwl, bwu+1 ) and use for
           off_diagonal block of reduced system.
         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} {{GU}_i} = {DL_i}$ .
         Transpose transmitted triangular matrix  $DL_i$
         Calculate the "spike" fillin, ${U_i}^T {{GL}_i}^T = {DU_i}^T$
         Copy D block into AF storage for solve.
         Calculate the update block for previous proc, E_i = GL_i{GU_i}
         Zero out space in case result is smaller than storage block
         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.
           Copy matrix HU_i (the last bwl rows of GU_i) to AFL storage
             as per requirements of BLAS routine DTRMM.
             Since we have GU_i stored,
             transpose HU_i to HU_i^T.
           Copy matrix HL_i (the last bwu rows of GL_i^T) to AFU store
             as per requirements of BLAS routine DTRMM.
             Since we have GL_i^T stored,
             transpose HL_i^T to HL_i.
       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'...
       Factor diagonal block
       ****************************************************************
       Receive offdiagonal block from processor to right.
         If this is the first group of processors, the receive comes
         from a different processor than otherwise.
 
Display dynamic version Find AutoScroll Reload FontSize: - + Hide Comments Hide Blanks Frame FullScreen MailPrint

  
001        SUBROUTINE PSDBTRF( N , BWL , BWU , A , JA , DESCA , AF , LAF , WORK ,
002       $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  *     March 12 , 2002
008  
009  *     .. Scalar Arguments ..
010        INTEGER BWL , BWU , INFO , JA , LAF , LWORK , N
011  *     ..
012  *     .. Array Arguments ..
013        INTEGER DESCA( * )
014        REAL A( * ) , AF( * ) , WORK( * )
015  *     ..
016  
017  *     Purpose
018  *     === ====
019  
020  *     PSDBTRF computes a LU factorization
021  *     of an N - by - N real banded
022  *     diagonally dominant - like distributed matrix
023  *     with bandwidth BWL , BWU : 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 PSDBTRS 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 banded upper triangular matrix and L is banded
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  *     BWL(global input) INTEGER
047  *     Number of subdiagonals. 0 <= BWL <= N - 1
048  
049  *     BWU(global input) INTEGER
050  *     Number of superdiagonals. 0 <= BWU <= N - 1
051  
052  *     A(local input / local output) REAL pointer into
053  *     local memory to an array with first dimension
054  *     LLD_A >=(bwl + bwu + 1)(stored in DESCA).
055  *     On entry , this array contains the local pieces of the
056  *     N - by - N unsymmetric banded distributed matrix
057  *     A(1 : N , JA : JA + N - 1) to be factored.
058  *     This local portion is stored in the packed banded format
059  *     used in LAPACK. Please see the Notes below and the
060  *     ScaLAPACK manual for more detail on the format of
061  *     distributed matrices.
062  *     On exit , this array contains information containing details
063  *     of the factorization.
064  *     Note that permutations are performed on the matrix , so that
065  *     the factors returned are different from those returned
066  *     by LAPACK.
067  
068  *     JA(global input) INTEGER
069  *     The index in the global array A that points to the start of
070  *     the matrix to be operated on(which may be either all of A
071  *     or a submatrix of A).
072  
073  *     DESCA(global and local input) INTEGER array of dimension DLEN.
074  *     if 1D type(DTYPE_A = 501) , DLEN >= 7 ;
075  *     if 2D type(DTYPE_A = 1) , DLEN >= 9 .
076  *     The array descriptor for the distributed matrix A.
077  *     Contains information of mapping of A to memory. Please
078  *     see NOTES below for full description and options.
079  
080  *     AF(local output) REAL array , dimension LAF.
081  *     Auxiliary Fillin Space.
082  *     Fillin is created during the factorization routine
083  *     PSDBTRF and this is stored in AF. If a linear system
084  *     is to be solved using PSDBTRS after the factorization
085  *     routine , AF *must not be altered* after the factorization.
086  
087  *     LAF(local input) INTEGER
088  *     Size of user - input Auxiliary Fillin space AF. Must be >=
089  *     NB*(bwl + bwu) + 6*max(bwl , bwu)*max(bwl , bwu)
090  *     If LAF is not large enough , an error code will be returned
091  *     and the minimum acceptable size will be returned in AF( 1 )
092  
093  *     WORK(local workspace / local output)
094  *     REAL temporary workspace. This space may
095  *     be overwritten in between calls to routines. WORK must be
096  *     the size given in LWORK.
097  *     On exit , WORK( 1 ) contains the minimal LWORK.
098  
099  *     LWORK(local input or global input) INTEGER
100  *     Size of user - input workspace WORK.
101  *     If LWORK is too small , the minimal acceptable size will be
102  *     returned in WORK(1) and an error code is returned. LWORK >=
103  *     max(bwl , bwu)*max(bwl , bwu)
104  
105  *     INFO(global output) INTEGER
106  *     = 0 : successful exit
107  *     < 0 : If the i - th argument is an array and the j - entry had
108  *     an illegal value , then INFO = - (i*100 + j) , if the i - th
109  *     argument is a scalar and had an illegal value , then
110  *     INFO = - i.
111  *     > 0 : If INFO = K <= NPROCS , the submatrix stored on processor
112  *     INFO and factored locally was not
113  *     diagonally dominant - like , and
114  *     the factorization was not completed.
115  *     If INFO = K > NPROCS , the submatrix stored on processor
116  *     INFO - NPROCS representing interactions with other
117  *     processors was not
118  *     stably factorable wo / interchanges ,
119  *     and the factorization was not completed.
120  
121  *     === ==================================================================
122  
123  *     Restrictions
124  *     === =========
125  
126  *     The following are restrictions on the input parameters. Some of these
127  *     are temporary and will be removed in future releases , while others
128  *     may reflect fundamental technical limitations.
129  
130  *     Non - cyclic restriction : VERY IMPORTANT !
131  *     P*NB >= mod(JA - 1 , NB) + N.
132  *     The mapping for matrices must be blocked , reflecting the nature
133  *     of the divide and conquer algorithm as a task - parallel algorithm.
134  *     This formula in words is : no processor may have more than one
135  *     chunk of the matrix.
136  
137  *     Blocksize cannot be too small :
138  *     If the matrix spans more than one processor , the following
139  *     restriction on NB , the size of each block on each processor ,
140  *     must hold :
141  *     NB >= 2*MAX(BWL , BWU)
142  *     The bulk of parallel computation is done on the matrix of size
143  *     O(NB) on each processor. If this is too small , divide and conquer
144  *     is a poor choice of algorithm.
145  
146  *     Submatrix reference :
147  *     JA = IB
148  *     Alignment restriction that prevents unnecessary communication.
149  
150  *     === ==================================================================
151  
152  *     Notes
153  *     === ==
154  
155  *     If the factorization routine and the solve routine are to be called
156  *     separately(to solve various sets of righthand sides using the same
157  *     coefficient matrix) , the auxiliary space AF *must not be altered*
158  *     between calls to the factorization routine and the solve routine.
159  
160  *     The best algorithm for solving banded and tridiagonal linear systems
161  *     depends on a variety of parameters , especially the bandwidth.
162  *     Currently , only algorithms designed for the case N / P >> bw are
163  *     implemented. These go by many names , including Divide and Conquer ,
164  *     Partitioning , domain decomposition - type , etc.
165  
166  *     Algorithm description : Divide and Conquer
167  
168  *     The Divide and Conqer algorithm assumes the matrix is narrowly
169  *     banded compared with the number of equations. In this situation ,
170  *     it is best to distribute the input matrix A one - dimensionally ,
171  *     with columns atomic and rows divided amongst the processes.
172  *     The basic algorithm divides the banded matrix up into
173  *     P pieces with one stored on each processor ,
174  *     and then proceeds in 2 phases for the factorization or 3 for the
175  *     solution of a linear system.
176  *     1) Local Phase :
177  *     The individual pieces are factored independently and in
178  *     parallel. These factors are applied to the matrix creating
179  *     fillin , which is stored in a non - inspectable way in auxiliary
180  *     space AF. Mathematically , this is equivalent to reordering
181  *     the matrix A as P A P^T and then factoring the principal
182  *     leading submatrix of size equal to the sum of the sizes of
183  *     the matrices factored on each processor. The factors of
184  *     these submatrices overwrite the corresponding parts of A
185  *     in memory.
186  *     2) Reduced System Phase :
187  *     A small(max(bwl , bwu)* (P - 1)) system is formed representing
188  *     interaction of the larger blocks , and is stored(as are its
189  *     factors) in the space AF. A parallel Block Cyclic Reduction
190  *     algorithm is used. For a linear system , a parallel front solve
191  *     followed by an analagous backsolve , both using the structure
192  *     of the factored matrix , are performed.
193  *     3) Backsubsitution Phase :
194  *     For a linear system , a local backsubstitution is performed on
195  *     each processor in parallel.
196  
197  *     Descriptors
198  *     === ========
199  
200  *     Descriptors now have *types* and differ from ScaLAPACK 1.0.
201  
202  *     Note : banded codes can use either the old two dimensional
203  *     or new one - dimensional descriptors , though the processor grid in
204  *     both cases *must be one - dimensional*. We describe both types below.
205  
206  *     Each global data object is described by an associated description
207  *     vector. This vector stores the information required to establish
208  *     the mapping between an object element and its corresponding process
209  *     and memory location.
210  
211  *     Let A be a generic term for any 2D block cyclicly distributed array.
212  *     Such a global array has an associated description vector DESCA.
213  *     In the following comments , the character _ should be read as
214  *     "of the global array".
215  
216  *     NOTATION STORED IN EXPLANATION
217  *     --- ------------ -------------- --------------------------------------
218  *     DTYPE_A(global) DESCA( DTYPE_ )The descriptor type. In this case ,
219  *     DTYPE_A = 1.
220  *     CTXT_A(global) DESCA( CTXT_ ) The BLACS context handle , indicating
221  *     the BLACS process grid A is distribu -
222  *     ted over. The context itself is glo -
223  *     bal , but the handle(the integer
224  *     value) may vary.
225  *     M_A(global) DESCA( M_ ) The number of rows in the global
226  *     array A.
227  *     N_A(global) DESCA( N_ ) The number of columns in the global
228  *     array A.
229  *     MB_A(global) DESCA( MB_ ) The blocking factor used to distribute
230  *     the rows of the array.
231  *     NB_A(global) DESCA( NB_ ) The blocking factor used to distribute
232  *     the columns of the array.
233  *     RSRC_A(global) DESCA( RSRC_ ) The process row over which the first
234  *     row of the array A is distributed.
235  *     CSRC_A(global) DESCA( CSRC_ ) The process column over which the
236  *     first column of the array A is
237  *     distributed.
238  *     LLD_A(local) DESCA( LLD_ ) The leading dimension of the local
239  *     array. LLD_A >= MAX(1 , LOCr(M_A)).
240  
241  *     Let K be the number of rows or columns of a distributed matrix ,
242  *     and assume that its process grid has dimension p x q.
243  *     LOCr( K ) denotes the number of elements of K that a process
244  *     would receive if K were distributed over the p processes of its
245  *     process column.
246  *     Similarly , LOCc( K ) denotes the number of elements of K that a
247  *     process would receive if K were distributed over the q processes of
248  *     its process row.
249  *     The values of LOCr() and LOCc() may be determined via a call to the
250  *     ScaLAPACK tool function , NUMROC :
251  *     LOCr( M ) = NUMROC( M , MB_A , MYROW , RSRC_A , NPROW ) ,
252  *     LOCc( N ) = NUMROC( N , NB_A , MYCOL , CSRC_A , NPCOL ).
253  *     An upper bound for these quantities may be computed by :
254  *     LOCr( M ) <= ceil( ceil(M / MB_A) / NPROW )*MB_A
255  *     LOCc( N ) <= ceil( ceil(N / NB_A) / NPCOL )*NB_A
256  
257  *     One - dimensional descriptors :
258  
259  *     One - dimensional descriptors are a new addition to ScaLAPACK since
260  *     version 1.0. They simplify and shorten the descriptor for 1D
261  *     arrays.
262  
263  *     Since ScaLAPACK supports two - dimensional arrays as the fundamental
264  *     object , we allow 1D arrays to be distributed either over the
265  *     first dimension of the array(as if the grid were P - by - 1) or the
266  *     2nd dimension(as if the grid were 1 - by - P). This choice is
267  *     indicated by the descriptor type(501 or 502)
268  *     as described below.
269  
270  *     IMPORTANT NOTE : the actual BLACS grid represented by the
271  *     CTXT entry in the descriptor may be *either* P - by - 1 or 1 - by - P
272  *     irrespective of which one - dimensional descriptor type
273  *     (501 or 502) is input.
274  *     This routine will interpret the grid properly either way.
275  *     ScaLAPACK routines *do not support intercontext operations* so that
276  *     the grid passed to a single ScaLAPACK routine *must be the same*
277  *     for all array descriptors passed to that routine.
278  
279  *     NOTE : In all cases where 1D descriptors are used , 2D descriptors
280  *     may also be used , since a one - dimensional array is a special case
281  *     of a two - dimensional array with one dimension of size unity.
282  *     The two - dimensional array used in this case *must* be of the
283  *     proper orientation :
284  *     If the appropriate one - dimensional descriptor is DTYPEA = 501
285  *     (1 by P type) , then the two dimensional descriptor must
286  *     have a CTXT value that refers to a 1 by P BLACS grid ;
287  *     If the appropriate one - dimensional descriptor is DTYPEA = 502
288  *     (P by 1 type) , then the two dimensional descriptor must
289  *     have a CTXT value that refers to a P by 1 BLACS grid.
290  
291  *     Summary of allowed descriptors , types , and BLACS grids :
292  *     DTYPE 501 502 1 1
293  *     BLACS grid 1xP or Px1 1xP or Px1 1xP Px1
294  *     --- --------------------------------------------------
295  *     A               OK NO OK NO
296  *     B               NO OK NO OK
297  
298  *     Let A be a generic term for any 1D block cyclicly distributed array.
299  *     Such a global array has an associated description vector DESCA.
300  *     In the following comments , the character _ should be read as
301  *     "of the global array".
302  
303  *     NOTATION STORED IN EXPLANATION
304  *     --- ------------ ---------- ------------------------------------------
305  *     DTYPE_A(global) DESCA( 1 ) The descriptor type. For 1D grids ,
306  *     TYPE_A = 501 : 1 - by - P grid.
307  *     TYPE_A = 502 : P - by - 1 grid.
308  *     CTXT_A(global) DESCA( 2 ) The BLACS context handle , indicating
309  *     the BLACS process grid A is distribu -
310  *     ted over. The context itself is glo -
311  *     bal , but the handle(the integer
312  *     value) may vary.
313  *     N_A(global) DESCA( 3 ) The size of the array dimension being
314  *     distributed.
315  *     NB_A(global) DESCA( 4 ) The blocking factor used to distribute
316  *     the distributed dimension of the array.
317  *     SRC_A(global) DESCA( 5 ) The process row or column over which the
318  *     first row or column of the array
319  *     is distributed.
320  *     LLD_A(local) DESCA( 6 ) The leading dimension of the local array
321  *     storing the local blocks of the distri -
322  *     buted array A. Minimum value of LLD_A
323  *     depends on TYPE_A.
324  *     TYPE_A = 501 : LLD_A >=
325  *     size of undistributed dimension , 1.
326  *     TYPE_A = 502 : LLD_A >= NB_A , 1.
327  *     Reserved DESCA( 7 ) Reserved for future use.
328  
329  *     === ==================================================================
330  
331  *     Code Developer : Andrew J. Cleary , University of Tennessee.
332  *     Current address : Lawrence Livermore National Labs.
333  *     Last modified by : Peter Arbenz , Institute of Scientific Computing ,
334  *     ETH , Zurich.
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        IF( LEVEL_DIST.EQ.1 ) THEN
353            COMM_PROC = MYCOL + 1
354  
355  *         Move block into place that it will be expected to be for
356  *         calcs.
357  
358            CALL SLACPY( 'N' , MAX_BW , MAX_BW , AF( ODD_SIZE*BWU + 1 ) , MAX_BW ,
359       $    AF( WORK_U + ODD_SIZE*BWL + 2*MBW2 + 1 ) , MAX_BW )
360  
361            CALL SLACPY( 'N' , MAX_BW , MAX_BW , AF( WORK_U + ODD_SIZE*BWL + 1 ) ,
362       $    MAX_BW , AF( ODD_SIZE*BWU + 2*MBW2 + 1 ) , MAX_BW )
363  
364        ELSE
365            COMM_PROC = MYCOL + LEVEL_DIST / 2
366        END IF
367  
368        IF( MYCOL / LEVEL_DIST.LE.( NPCOL - 1 ) / LEVEL_DIST - 2 ) THEN
369  
370            CALL SGERV2D( ICTXT , MAX_BW , MAX_BW , AF( ODD_SIZE*BWU + 1 ) ,
371       $    MAX_BW , 0 , COMM_PROC )
372  
373            CALL SGERV2D( ICTXT , MAX_BW , MAX_BW ,
374       $    AF( WORK_U + ODD_SIZE*BWL + 1 ) , MAX_BW , 0 ,
375       $    COMM_PROC )
376  
377            IF( INFO.EQ.0 ) THEN
378  
379  *             Modify upper off_diagonal block with diagonal block
380  
381                CALL STBTRS( 'L' , 'N' , 'U' , BWU , MIN( BWL , BWU - 1 ) , BWU ,
382       $        AF( ODD_SIZE*BWU + MBW2 + 1 + ( MAX_BW + 1 )*( MAX_BW -
383       $        BWU ) ) , MAX_BW + 1 , AF( WORK_U + ODD_SIZE*BWL + 1 +
384       $        MAX_BW - BWU ) , MAX_BW , INFO )
385  
386  *             Modify lower off_diagonal block with diagonal block
387  
388                CALL STBTRS( 'U' , 'T' , 'N' , BWL , MIN( BWU , BWL - 1 ) , BWL ,
389       $        AF( ODD_SIZE*BWU + MBW2 + 1 - MIN( BWU ,
390       $        BWL - 1 ) + ( MAX_BW + 1 )*( MAX_BW - BWL ) ) ,
391       $        MAX_BW + 1 , AF( ODD_SIZE*BWU + 1 + MAX_BW - BWL ) ,
392       $        MAX_BW , INFO )
393  
394            END IF
395  *         End of "if( info.eq.0 ) then"
396  
397  *         Calculate contribution from this block to next diagonal block
398  
399            CALL SGEMM( 'T' , 'N' , MAX_BW , MAX_BW , MAX_BW , - ONE ,
400       $    AF(( ODD_SIZE )*BWU + 1 ) , MAX_BW ,
401       $    AF( WORK_U + ( ODD_SIZE )*BWL + 1 ) , MAX_BW , ZERO ,
402       $    WORK( 1 ) , MAX_BW )
403  
404  *         Send contribution to diagonal block's owning processor.
405  
406            CALL SGESD2D( ICTXT , MAX_BW , MAX_BW , WORK( 1 ) , MAX_BW , 0 ,
407       $    MYCOL + LEVEL_DIST )
408  
409        END IF
410  *     End of "if( mycol/level_dist .le.(npcol-1)/level_dist-2 )..."
411  
412  *     ****************************************************************
413  *     Receive off_diagonal block from left and use to finish with this
414  *     processor.
415  
416        IF(( MYCOL / LEVEL_DIST.GT.0 ) .AND.
417       $( MYCOL / LEVEL_DIST.LE.( NPCOL - 1 ) / LEVEL_DIST - 1 ) ) THEN
418  
419        IF( LEVEL_DIST.GT.1 ) THEN
420  
421  *         Receive offdiagonal block(s) from proc level_dist / 2 to the
422  *         left
423  
424            CALL SGERV2D( ICTXT , MAX_BW , MAX_BW ,
425       $    AF( WORK_U + ODD_SIZE*BWL + 2*MBW2 + 1 ) , MAX_BW , 0 ,
426       $    MYCOL - LEVEL_DIST / 2 )
427  
428  *         Receive offdiagonal block(s) from proc level_dist / 2 to the
429  *         left
430  
431            CALL SGERV2D( ICTXT , MAX_BW , MAX_BW ,
432       $    AF( ODD_SIZE*BWU + 2*MBW2 + 1 ) , MAX_BW , 0 ,
433       $    MYCOL - LEVEL_DIST / 2 )
434  
435        END IF
436  
437        IF( INFO.EQ.0 ) THEN
438  
439  *         Use diagonal block(s) to modify this offdiagonal block
440  
441  *         Since STBTRS has no "left-right" option , we must transpose
442  
443            CALL SLATCPY( 'N' , MAX_BW , MAX_BW ,
444       $    AF( WORK_U + ODD_SIZE*BWL + 2*MBW2 + 1 ) , MAX_BW ,
445       $    WORK( 1 ) , MAX_BW )
446  
447            CALL STBTRS( 'L' , 'N' , 'U' , MAX_BW , MIN( BWL , MAX_BW - 1 ) ,
448       $    BWL , AF( ODD_SIZE*BWU + MBW2 + 1 ) , MAX_BW + 1 ,
449       $    WORK( 1 + MAX_BW*( MAX_BW - BWL ) ) , MAX_BW , INFO )
450  
451  *         Transpose back
452  
453            CALL SLATCPY( 'N' , MAX_BW , MAX_BW , WORK( 1 ) , MAX_BW ,
454       $    AF( WORK_U + ODD_SIZE*BWL + 2*MBW2 + 1 ) , MAX_BW )
455  
456  *         Since STBTRS has no "left-right" option , we must transpose
457  
458            CALL SLATCPY( 'N' , MAX_BW , MAX_BW ,
459       $    AF( ODD_SIZE*BWU + 2*MBW2 + 1 ) , MAX_BW ,
460       $    WORK( 1 ) , MAX_BW )
461  
462            CALL STBTRS( 'U' , 'T' , 'N' , MAX_BW , MIN( BWU , MAX_BW - 1 ) ,
463       $    BWU , AF( ODD_SIZE*BWU + MBW2 + 1 - MIN( BWU ,
464       $    MAX_BW - 1 ) ) , MAX_BW + 1 ,
465       $    WORK( 1 + MAX_BW*( MAX_BW - BWU ) ) , MAX_BW , INFO )
466  
467  *         Transpose back
468  
469            CALL SLATCPY( 'N' , MAX_BW , MAX_BW , WORK( 1 ) , MAX_BW ,
470       $    AF( ODD_SIZE*BWU + 2*MBW2 + 1 ) , MAX_BW )
471  
472        END IF
473  *     End of "if( info.eq.0 ) then"
474  
475  *     Use offdiag block(s) to calculate modification to diag block
476  *     of processor to the left
477  
478        CALL SGEMM( 'N' , 'T' , MAX_BW , MAX_BW , MAX_BW , - ONE ,
479       $AF(( ODD_SIZE )*BWU + 2*MBW2 + 1 ) , MAX_BW ,
480       $AF( WORK_U + ( ODD_SIZE )*BWL + 2*MBW2 + 1 ) , MAX_BW ,
481       $ZERO , WORK( 1 ) , MAX_BW )
482  
483  *     Send contribution to diagonal block's owning processor.
484  
485        CALL SGESD2D( ICTXT , MAX_BW , MAX_BW , WORK( 1 ) , MAX_BW , 0 ,
486       $MYCOL - LEVEL_DIST )
487  
488  *     *******************************************************
489  
490        IF( MYCOL / LEVEL_DIST.LE.( NPCOL - 1 ) / LEVEL_DIST - 2 ) THEN
491  
492  *         Decide which processor offdiagonal block(s) goes to
493  
494            IF(( MOD( MYCOL / ( 2*LEVEL_DIST ) , 2 ) ).EQ.0 ) THEN
495                COMM_PROC = MYCOL + LEVEL_DIST
496            ELSE
497                COMM_PROC = MYCOL - LEVEL_DIST
498            END IF
499  
500  *         Use offdiagonal blocks to calculate offdiag
501  *         block to send to neighboring processor. Depending
502  *         on circumstances , may need to transpose the matrix.
503  
504            CALL SGEMM( 'N' , 'N' , MAX_BW , MAX_BW , MAX_BW , - ONE ,
505       $    AF( WORK_U + ODD_SIZE*BWL + 2*MBW2 + 1 ) , MAX_BW ,
506       $    AF( ODD_SIZE*BWU + 1 ) , MAX_BW , ZERO , WORK( 1 ) ,
507       $    MAX_BW )
508  
509  *         Send contribution to offdiagonal block's owning processor.
510  
511            CALL SGESD2D( ICTXT , MAX_BW , MAX_BW , WORK( 1 ) , MAX_BW , 0 ,
512       $    COMM_PROC )
513  
514            CALL SGEMM( 'N' , 'N' , MAX_BW , MAX_BW , MAX_BW , - ONE ,
515       $    AF( ODD_SIZE*BWU + 2*MBW2 + 1 ) , MAX_BW ,
516       $    AF( WORK_U + ODD_SIZE*BWL + 1 ) , MAX_BW , ZERO ,
517       $    WORK( 1 ) , MAX_BW )
518  
519  *         Send contribution to offdiagonal block's owning processor.
520  
521            CALL SGESD2D( ICTXT , MAX_BW , MAX_BW , WORK( 1 ) , MAX_BW , 0 ,
522       $    COMM_PROC )
523  
524        END IF
525  
526        END IF
527  *     End of "if( mycol/level_dist.le.(npcol-1)/level_dist -1 )..."
528  
529    120 CONTINUE
530  
531    130 CONTINUE
532  
533  *     Free BLACS space used to hold standard - form grid.
534  
535        IF( ICTXT_SAVE.NE.ICTXT_NEW ) THEN
536            CALL BLACS_GRIDEXIT( ICTXT_NEW )
537        END IF
538  
539    140 CONTINUE
540  
541  *     Restore saved input parameters
542  
543        ICTXT = ICTXT_SAVE
544        NP = NP_SAVE
545  
546  *     Output minimum worksize
547  
548        WORK( 1 ) = WORK_SIZE_MIN
549  
550  *     Make INFO consistent across processors
551  
552        CALL IGAMX2D( ICTXT , 'A' , ' ' , 1 , 1 , INFO , 1 , INFO , INFO , - 1 , 0 ,
553       $0 )
554  
555        IF( MYCOL.EQ.0 ) THEN
556            CALL IGEBS2D( ICTXT , 'A' , ' ' , 1 , 1 , INFO , 1 )
557        ELSE
558            CALL IGEBR2D( ICTXT , 'A' , ' ' , 1 , 1 , INFO , 1 , 0 , 0 )
559        END IF
560  
561        RETURN
562  
563  *     End of PSDBTRF
564  
565        END