SRC\pslasrt.f

#lines: 254   size: 8 Kb   creation: 18/01/2006 23:36:04   last modification: 08/05/2008 18:38:06   attribute: ARCH    Find   Reload  
1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:
22:
23:
24:
25:
26:
27:
28:
29:
30:
31:
32:
33:
34:
35:
36:
37:
38:
39:
40:
41:
42:
43:
44:
45:
46:
47:
48:
49:
50:
51:
52:
53:
54:
55:
56:
57:
58:
59:
60:
61:
62:
63:
64:
65:
66:
67:
68:
69:
70:
71:
72:
73:
74:
75:
76:
77:
78:
79:
80:
81:
82:
83:
84:
85:
86:
87:
88:
89:
90:
91:
92:
93:
94:
95:
96:
97:
98:
99:
100:
101:
102:
103:
104:
105:
106:
107:
108:
109:
110:
111:
112:
113:
114:
115:
116:
117:
118:
119:
120:
121:
122:
123:
124:
125:
126:
127:
128:
129:
130:
131:
132:
133:
134:
135:
136:
137:
138:
139:
140:
141:
142:
143:
144:
145:
146:
147:
148:
149:
150:
151:
152:
153:
154:
155:
156:
157:
158:
159:
160:
161:
162:
163:
164:
165:
166:
167:
168:
169:
170:
171:
172:
173:
174:
175:
176:
177:
178:
179:
180:
181:
182:
183:
184:
185:
186:
187:
188:
189:
190:
191:
192:
193:
194:
195:
196:
197:
198:
199:
200:
201:
202:
203:
204:
205:
206:
207:
208:
209:
210:
211:
212:
213:
214:
215:
216:
217:
218:
219:
220:
221:
222:
223:
224:
225:
226:
227:
228:
229:
230:
231:
232:
233:
234:
235:
236:
237:
238:
239:
240:
241:
242:
243:
244:
245:
246:
247:
248:
249:
250:
251:
252:
253:
254:
      SUBROUTINE PSLASRT( ID, N, D, Q, IQ, JQ, DESCQ, WORK, LWORK, 
     $                    IWORK, LIWORK, INFO )
*
*  -- ScaLAPACK auxiliary routine (version 1.7) --
*     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
*     and University of California, Berkeley.
*     February 22, 2000
*
*     .. Scalar Arguments ..
      CHARACTER          ID
      INTEGER            INFO, IQ, JQ, LIWORK, LWORK, N
*     ..
*     .. Array Arguments ..
      INTEGER            DESCQ( * ), IWORK( * )
      REAL               D( * ), Q( * ), WORK( * )
*     ..
*
*  Purpose
*  =======
*
*  PSLASRT Sort the numbers in D in increasing order and the
*  corresponding vectors in Q.
*
*  Arguments
*  =========
*
*  ID      (global input) CHARACTER*1
*          = 'I': sort D in increasing order;
*          = 'D': sort D in decreasing order. (NOT IMPLEMENTED YET)
*
*  N       (global input) INTEGER
*          The number of columns to be operated on i.e the number of
*          columns of the distributed submatrix sub( Q ). N >= 0.
*
*  D       (global input/output) REAL array, dimmension (N)
*          On exit, the number in D are sorted in increasing order.
*
*  Q       (local input) REAL pointer into the local memory
*          to an array of dimension (LLD_Q, LOCc(JQ+N-1) ). This array
*          contains the local pieces of the distributed matrix sub( A )
*          to be copied from.
*
*  IQ      (global input) INTEGER
*          The row index in the global array A indicating the first
*          row of sub( Q ).
*
*  JQ      (global input) INTEGER
*          The column index in the global array A indicating the
*          first column of sub( Q ).
*
*  DESCQ   (global and local input) INTEGER array of dimension DLEN_.
*          The array descriptor for the distributed matrix A.
*
*  WORK    (local workspace/local output) REAL array,
*                                                  dimension (LWORK)
*  LWORK   (local or global input) INTEGER
*          The dimension of the array WORK.
*          LWORK = MAX( N, NP * ( NB + NQ ))
*          where
*          NP = NUMROC( N, NB, MYROW, IAROW, NPROW ),
*          NQ = NUMROC( N, NB, MYCOL, DESCQ( CSRC_ ), NPCOL )
*
*  IWORK   (local workspace/local output) INTEGER array,
*                                                  dimension (LIWORK)
*
*  LIWORK (local or global input) INTEGER
*          The dimension of the array IWORK.
*          LIWORK = N + 2*NB + 2*NPCOL
*
*  INFO    (global output) INTEGER
*          = 0:  successful exit
*          < 0:  If the i-th argument is an array and the j-entry had
*                an illegal value, then INFO = -(i*100+j), if the i-th
*                argument is a scalar and had an illegal value, then
*                INFO = -i.
*
*  =====================================================================
*
*     .. Parameters ..
      INTEGER            BLOCK_CYCLIC_2D, DLEN_, DTYPE_, CTXT_, M_, N_,
     $                   MB_, NB_, RSRC_, CSRC_, LLD_
      PARAMETER          ( BLOCK_CYCLIC_2D = 1, DLEN_ = 9, DTYPE_ = 1,
     $                   CTXT_ = 2, M_ = 3, N_ = 4, MB_ = 5, NB_ = 6,
     $                   RSRC_ = 7, CSRC_ = 8, LLD_ = 9 )
*     ..
*     .. Local Scalars ..
      INTEGER            CL, COL, DUMMY, I, ICTXT, IID, IIQ, INDCOL,
     $                   INDX, INDXC, INDXG, IPQ, IPQ2, IPW, IPWORK, J,
     $                   JJQ, K, L, LDQ, LEND, LIWMIN, LWMIN, MYCOL,
     $                   MYROW, NB, ND, NP, NPCOL, NPROW, NQ, PSQ, QCOL,
     $                   QTOT, SBUF
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      INTEGER            INDXG2L, INDXG2P, NUMROC
      EXTERNAL           INDXG2L, INDXG2P, LSAME, NUMROC
*     ..
*     .. External Subroutines ..
      EXTERNAL           BLACS_GRIDINFO, CHK1MAT, PXERBLA, SCOPY,
     $                   SGERV2D, SGESD2D, SLACPY, SLAPST
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          MAX, MIN, MOD
*     ..
*     .. Executable Statements ..
*
*       This is just to keep ftnchek and toolpack/1 happy
      IF( BLOCK_CYCLIC_2D*CSRC_*CTXT_*DLEN_*DTYPE_*LLD_*MB_*M_*NB_*N_*
     $    RSRC_.LT.0 )RETURN
*
      IF( N.EQ.0 )
     $   RETURN
*
      CALL BLACS_GRIDINFO( DESCQ( CTXT_ ), NPROW, NPCOL, MYROW, MYCOL )
*
*     Test the input parameters
*
      INFO = 0
      IF( NPROW.EQ.-1 ) THEN
         INFO = -( 600+CTXT_ )
      ELSE
         CALL CHK1MAT( N, 1, N, 1, IQ, JQ, DESCQ, 6, INFO )
         IF( INFO.EQ.0 ) THEN
            NB = DESCQ( NB_ )
            LDQ = DESCQ( LLD_ )
            NP = NUMROC( N, NB, MYROW, DESCQ( RSRC_ ), NPROW )
            NQ = NUMROC( N, NB, MYCOL, DESCQ( CSRC_ ), NPCOL )
            LWMIN = MAX( N, NP*( NB+NQ ) )
            LIWMIN = N + 2*( NB+NPCOL )
            IF( .NOT.LSAME( ID, 'I' ) ) THEN
               INFO = -1
            ELSE IF( N.LT.0 ) THEN
               INFO = -2
            ELSE IF( LWORK.LT.LWMIN ) THEN
               INFO = -9
            ELSE IF( LIWORK.LT.LIWMIN ) THEN
               INFO = -11
            END IF
         END IF
      END IF
*
      IF( INFO.NE.0 ) THEN
         CALL PXERBLA( ICTXT, 'PSLASRT', -INFO )
         RETURN
      END IF
*
*     Set Pointers
*
      INDXC = 1
      INDX = INDXC + N
      INDXG = INDX
      INDCOL = INDXG + NB
      QTOT = INDCOL + NB
      PSQ = QTOT + NPCOL
*
      IID = 1
      IPQ2 = 1
      IPW = IPQ2 + NP*NQ
*
      DUMMY = 0
      IIQ = INDXG2L( IQ, NB, DUMMY, DUMMY, NPROW )
*
*     Sort the eigenvalues in D
*
      CALL SLAPST( 'I', N, D, IWORK( INDX ), INFO )
*
      DO 10 L = 0, N - 1
         WORK( IID+L ) = D( IWORK( INDX+L ) )
         IWORK( INDXC-1+IWORK( INDX+L ) ) = IID + L
   10 CONTINUE
      CALL SCOPY( N, WORK, 1, D, 1 )
*
      ND = 0
   20 CONTINUE
      IF( ND.LT.N ) THEN
         LEND = MIN( NB, N-ND )
         J = JQ + ND
         QCOL = INDXG2P( J, NB, DUMMY, DESCQ( CSRC_ ), NPCOL )
         K = 0
         DO 30 L = 0, LEND - 1
            I = JQ - 1 + IWORK( INDXC+ND+L )
            CL = INDXG2P( I, NB, DUMMY, DESCQ( CSRC_ ), NPCOL )
            IWORK( INDCOL+L ) = CL
            IF( MYCOL.EQ.CL ) THEN
               IWORK( INDXG+K ) = IWORK( INDXC+ND+L )
               K = K + 1
            END IF
   30    CONTINUE
*
         IF( MYCOL.EQ.QCOL ) THEN
            DO 40 CL = 0, NPCOL - 1
               IWORK( QTOT+CL ) = 0
   40       CONTINUE
            DO 50 L = 0, LEND - 1
               IWORK( QTOT+IWORK( INDCOL+L ) ) = IWORK( QTOT+
     $            IWORK( INDCOL+L ) ) + 1
   50       CONTINUE
            IWORK( PSQ ) = 1
            DO 60 CL = 1, NPCOL - 1
               IWORK( PSQ+CL ) = IWORK( PSQ+CL-1 ) + IWORK( QTOT+CL-1 )
   60       CONTINUE
            DO 70 L = 0, LEND - 1
               CL = IWORK( INDCOL+L )
               I = JQ + ND + L
               JJQ = INDXG2L( I, NB, DUMMY, DUMMY, NPCOL )
               IPQ = IIQ + ( JJQ-1 )*LDQ
               IPWORK = IPW + ( IWORK( PSQ+CL )-1 )*NP
               CALL SCOPY( NP, Q( IPQ ), 1, WORK( IPWORK ), 1 )
               IWORK( PSQ+CL ) = IWORK( PSQ+CL ) + 1
   70       CONTINUE
            IWORK( PSQ ) = 1
            DO 80 CL = 1, NPCOL - 1
               IWORK( PSQ+CL ) = IWORK( PSQ+CL-1 ) + IWORK( QTOT+CL-1 )
   80       CONTINUE
            DO 90 L = 0, K - 1
               I = IWORK( INDXG+L )
               JJQ = INDXG2L( I, NB, DUMMY, DUMMY, NPCOL )
               IPQ = IPQ2 + ( JJQ-1 )*NP
               IPWORK = IPW + ( IWORK( PSQ+MYCOL )-1 )*NP
               CALL SCOPY( NP, WORK( IPWORK ), 1, WORK( IPQ ), 1 )
               IWORK( PSQ+MYCOL ) = IWORK( PSQ+MYCOL ) + 1
   90       CONTINUE
            DO 100 CL = 1, NPCOL - 1
               COL = MOD( MYCOL+CL, NPCOL )
               SBUF = IWORK( QTOT+COL )
               IF( SBUF.NE.0 ) THEN
                  IPWORK = IPW + ( IWORK( PSQ+COL )-1 )*NP
                  CALL SGESD2D( DESCQ( CTXT_ ), NP, SBUF,
     $                          WORK( IPWORK ), NP, MYROW, COL )
               END IF
  100       CONTINUE
*
         ELSE
*
            IF( K.NE.0 ) THEN
               CALL SGERV2D( DESCQ( CTXT_ ), NP, K, WORK( IPW ), NP,
     $                       MYROW, QCOL )
               DO 110 L = 0, K - 1
                  I = JQ - 1 + IWORK( INDXG+L )
                  JJQ = INDXG2L( I, NB, DUMMY, DUMMY, NPCOL )
                  IPQ = 1 + ( JJQ-1 )*NP
                  IPWORK = IPW + L*NP
                  CALL SCOPY( NP, WORK( IPWORK ), 1, WORK( IPQ ), 1 )
  110          CONTINUE
            END IF
         END IF
         ND = ND + NB
         GO TO 20
      END IF
      CALL SLACPY( 'Full', NP, NQ, WORK, NP, Q( IIQ ), LDQ )
*
*     End of PSLASRT
*
      END