Context Navigation

source: trunk/SOURCES/band.f @ 25

Last change on this file since 25 was 4, checked in by dumas, 10 years ago
initial import GRISLI trunk
File size: 48.7 KB

Line
1
2	SUBROUTINE SGBSV( N, KL, KU, NRHS, AB, LDAB, IPIV, B, LDB, INFO )
3	*
4	* -- LAPACK driver routine (version 2.0) --
5	* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
6	* Courant Institute, Argonne National Lab, and Rice University
7	* March 31, 1993
8	*
9	* .. Scalar Arguments ..
10	INTEGER INFO, KL, KU, LDAB, LDB, N, NRHS
11	* ..
12	* .. Array Arguments ..
13	INTEGER IPIV( * )
14	REAL AB( LDAB, * ), B( LDB, * )
15	* ..
16	*
17	* Purpose
18	* =======
19	*
20	* SGBSV computes the solution to a real system of linear equations
21	* A * X = B, where A is a band matrix of order N with KL subdiagonals
22	* and KU superdiagonals, and X and B are N-by-NRHS matrices.
23	*
24	* The LU decomposition with partial pivoting and row interchanges is
25	* used to factor A as A = L * U, where L is a product of permutation
26	* and unit lower triangular matrices with KL subdiagonals, and U is
27	* upper triangular with KL+KU superdiagonals. The factored form of A
28	* is then used to solve the system of equations A * X = B.
29	*
30	* Arguments
31	* =========
32	*
33	* N (input) INTEGER
34	* The number of linear equations, i.e., the order of the
35	* matrix A. N >= 0.
36	*
37	* KL (input) INTEGER
38	* The number of subdiagonals within the band of A. KL >= 0.
39	*
40	* KU (input) INTEGER
41	* The number of superdiagonals within the band of A. KU >= 0.
42	*
43	* NRHS (input) INTEGER
44	* The number of right hand sides, i.e., the number of columns
45	* of the matrix B. NRHS >= 0.
46	*
47	* AB (input/output) REAL array, dimension (LDAB,N)
48	* On entry, the matrix A in band storage, in rows KL+1 to
49	* 2*KL+KU+1; rows 1 to KL of the array need not be set.
50	* The j-th column of A is stored in the j-th column of the
51	* array AB as follows:
52	* AB(KL+KU+1+i-j,j) = A(i,j) for max(1,j-KU)<=i<=min(N,j+KL)
53	* On exit, details of the factorization: U is stored as an
54	* upper triangular band matrix with KL+KU superdiagonals in
55	* rows 1 to KL+KU+1, and the multipliers used during the
56	* factorization are stored in rows KL+KU+2 to 2*KL+KU+1.
57	* See below for further details.
58	*
59	* LDAB (input) INTEGER
60	* The leading dimension of the array AB. LDAB >= 2*KL+KU+1.
61	*
62	* IPIV (output) INTEGER array, dimension (N)
63	* The pivot indices that define the permutation matrix P;
64	* row i of the matrix was interchanged with row IPIV(i).
65	*
66	* B (input/output) REAL array, dimension (LDB,NRHS)
67	* On entry, the N-by-NRHS right hand side matrix B.
68	* On exit, if INFO = 0, the N-by-NRHS solution matrix X.
69	*
70	* LDB (input) INTEGER
71	* The leading dimension of the array B. LDB >= max(1,N).
72	*
73	* INFO (output) INTEGER
74	* = 0: successful exit
75	* < 0: if INFO = -i, the i-th argument had an illegal value
76	* > 0: if INFO = i, U(i,i) is exactly zero. The factorization
77	* has been completed, but the factor U is exactly
78	* singular, and the solution has not been computed.
79	*
80	* Further Details
81	* ===============
82	*
83	* The band storage scheme is illustrated by the following example, when
84	* M = N = 6, KL = 2, KU = 1:
85	*
86	* On entry: On exit:
87	*
88	* * * * + + + * * * u14 u25 u36
89	* * * + + + + * * u13 u24 u35 u46
90	* * a12 a23 a34 a45 a56 * u12 u23 u34 u45 u56
91	* a11 a22 a33 a44 a55 a66 u11 u22 u33 u44 u55 u66
92	* a21 a32 a43 a54 a65 * m21 m32 m43 m54 m65 *
93	* a31 a42 a53 a64 * * m31 m42 m53 m64 * *
94	*
95	* Array elements marked * are not used by the routine; elements marked
96	* + need not be set on entry, but are required by the routine to store
97	* elements of U because of fill-in resulting from the row interchanges.
98	*
99	* =====================================================================
100	*
101	* .. External Subroutines ..
102	EXTERNAL SGBTRF, SGBTRS, XERBLA
103	* ..
104	* .. Intrinsic Functions ..
105	INTRINSIC MAX
106	* ..
107	* .. Executable Statements ..
108	*
109	* Test the input parameters.
110	*
111	c write(6,*) 'dans sgbsv'
112	c write(6,*) 'n',n,' kl',kl,' ku',ku,' nrhs',nrhs
113	c write(6,*) 'ldab',ldab,' ldb',ldb
114	INFO = 0
115	IF( N.LT.0 ) THEN
116	INFO = -1
117	ELSE IF( KL.LT.0 ) THEN
118	INFO = -2
119	ELSE IF( KU.LT.0 ) THEN
120	INFO = -3
121	ELSE IF( NRHS.LT.0 ) THEN
122	INFO = -4
123	ELSE IF( LDAB.LT.2*KL+KU+1 ) THEN
124	INFO = -6
125	ELSE IF( LDB.LT.MAX( N, 1 ) ) THEN
126	INFO = -9
127	END IF
128	IF( INFO.NE.0 ) THEN
129	CALL XERBLA( 'SGBSV ', -INFO )
130	RETURN
131	END IF
132	*
133	*
134	* Compute the LU factorization of the band matrix A.
135	*
136	CALL SGBTRF( N, N, KL, KU, AB, LDAB, IPIV, INFO )
137	IF( INFO.EQ.0 ) THEN
138	*
139	* Solve the system A*X = B, overwriting B with X.
140	*
141	CALL SGBTRS( 'No transpose', N, KL, KU, NRHS, AB, LDAB, IPIV,
142	$ B, LDB, INFO )
143	END IF
144	RETURN
145	*
146	* End of SGBSV
147	*
148	END
149	C =================================================================
150	SUBROUTINE SGBTF2( M, N, KL, KU, AB, LDAB, IPIV, INFO )
151	*
152	* -- LAPACK routine (version 2.0) --
153	* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
154	* Courant Institute, Argonne National Lab, and Rice University
155	* February 29, 1992
156	*
157	* .. Scalar Arguments ..
158	INTEGER INFO, KL, KU, LDAB, M, N
159	* ..
160	* .. Array Arguments ..
161	INTEGER IPIV( * )
162	REAL AB( LDAB, * )
163	* ..
164	*
165	* Purpose
166	* =======
167	*
168	* SGBTF2 computes an LU factorization of a real m-by-n band matrix A
169	* using partial pivoting with row interchanges.
170	*
171	* This is the unblocked version of the algorithm, calling Level 2 BLAS.
172	*
173	* Arguments
174	* =========
175	*
176	* M (input) INTEGER
177	* The number of rows of the matrix A. M >= 0.
178	*
179	* N (input) INTEGER
180	* The number of columns of the matrix A. N >= 0.
181	*
182	* KL (input) INTEGER
183	* The number of subdiagonals within the band of A. KL >= 0.
184	*
185	* KU (input) INTEGER
186	* The number of superdiagonals within the band of A. KU >= 0.
187	*
188	* AB (input/output) REAL array, dimension (LDAB,N)
189	* On entry, the matrix A in band storage, in rows KL+1 to
190	* 2*KL+KU+1; rows 1 to KL of the array need not be set.
191	* The j-th column of A is stored in the j-th column of the
192	* array AB as follows:
193	* AB(kl+ku+1+i-j,j) = A(i,j) for max(1,j-ku)<=i<=min(m,j+kl)
194	*
195	* On exit, details of the factorization: U is stored as an
196	* upper triangular band matrix with KL+KU superdiagonals in
197	* rows 1 to KL+KU+1, and the multipliers used during the
198	* factorization are stored in rows KL+KU+2 to 2*KL+KU+1.
199	* See below for further details.
200	*
201	* LDAB (input) INTEGER
202	* The leading dimension of the array AB. LDAB >= 2*KL+KU+1.
203	*
204	* IPIV (output) INTEGER array, dimension (min(M,N))
205	* The pivot indices; for 1 <= i <= min(M,N), row i of the
206	* matrix was interchanged with row IPIV(i).
207	*
208	* INFO (output) INTEGER
209	* = 0: successful exit
210	* < 0: if INFO = -i, the i-th argument had an illegal value
211	* > 0: if INFO = +i, U(i,i) is exactly zero. The factorization
212	* has been completed, but the factor U is exactly
213	* singular, and division by zero will occur if it is used
214	* to solve a system of equations.
215	*
216	* Further Details
217	* ===============
218	*
219	* The band storage scheme is illustrated by the following example, when
220	* M = N = 6, KL = 2, KU = 1:
221	*
222	* On entry: On exit:
223	*
224	* * * * + + + * * * u14 u25 u36
225	* * * + + + + * * u13 u24 u35 u46
226	* * a12 a23 a34 a45 a56 * u12 u23 u34 u45 u56
227	* a11 a22 a33 a44 a55 a66 u11 u22 u33 u44 u55 u66
228	* a21 a32 a43 a54 a65 * m21 m32 m43 m54 m65 *
229	* a31 a42 a53 a64 * * m31 m42 m53 m64 * *
230	*
231	* Array elements marked * are not used by the routine; elements marked
232	* + need not be set on entry, but are required by the routine to store
233	* elements of U, because of fill-in resulting from the row
234	* interchanges.
235	*
236	* =====================================================================
237	*
238	* .. Parameters ..
239	REAL ONE, ZERO
240	PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 )
241	* ..
242	* .. Local Scalars ..
243	INTEGER I, J, JP, JU, KM, KV
244	* ..
245	* .. External Functions ..
246	INTEGER ISAMAX
247	EXTERNAL ISAMAX
248	* ..
249	* .. External Subroutines ..
250	EXTERNAL SGER, SSCAL, SSWAP, XERBLA
251	* ..
252	* .. Intrinsic Functions ..
253	INTRINSIC MAX, MIN
254	* ..
255	* .. Executable Statements ..
256	*
257	* KV is the number of superdiagonals in the factor U, allowing for
258	* fill-in.
259	*
260	KV = KU + KL
261	*
262	* Test the input parameters.
263	*
264	INFO = 0
265	IF( M.LT.0 ) THEN
266	INFO = -1
267	ELSE IF( N.LT.0 ) THEN
268	INFO = -2
269	ELSE IF( KL.LT.0 ) THEN
270	INFO = -3
271	ELSE IF( KU.LT.0 ) THEN
272	INFO = -4
273	ELSE IF( LDAB.LT.KL+KV+1 ) THEN
274	INFO = -6
275	END IF
276	IF( INFO.NE.0 ) THEN
277	CALL XERBLA( 'SGBTF2', -INFO )
278	RETURN
279	END IF
280	*
281	* Quick return if possible
282	*
283	IF( M.EQ.0 .OR. N.EQ.0 )
284	$ RETURN
285	*
286	* Gaussian elimination with partial pivoting
287	*
288	* Set fill-in elements in columns KU+2 to KV to zero.
289	*
290	DO 20 J = KU + 2, MIN( KV, N )
291	DO 10 I = KV - J + 2, KL
292	AB( I, J ) = ZERO
293	10 CONTINUE
294	20 CONTINUE
295	*
296	* JU is the index of the last column affected by the current stage
297	* of the factorization.
298	*
299	JU = 1
300	*
301	DO 40 J = 1, MIN( M, N )
302	*
303	* Set fill-in elements in column J+KV to zero.
304	*
305	IF( J+KV.LE.N ) THEN
306	DO 30 I = 1, KL
307	AB( I, J+KV ) = ZERO
308	30 CONTINUE
309	END IF
310	*
311	* Find pivot and test for singularity. KM is the number of
312	* subdiagonal elements in the current column.
313	*
314	KM = MIN( KL, M-J )
315	JP = ISAMAX( KM+1, AB( KV+1, J ), 1 )
316	IPIV( J ) = JP + J - 1
317	IF( AB( KV+JP, J ).NE.ZERO ) THEN
318	JU = MAX( JU, MIN( J+KU+JP-1, N ) )
319	*
320	* Apply interchange to columns J to JU.
321	*
322	IF( JP.NE.1 )
323	$ CALL SSWAP( JU-J+1, AB( KV+JP, J ), LDAB-1,
324	$ AB( KV+1, J ), LDAB-1 )
325	*
326	IF( KM.GT.0 ) THEN
327	*
328	* Compute multipliers.
329	*
330	CALL SSCAL( KM, ONE / AB( KV+1, J ), AB( KV+2, J ), 1 )
331	*
332	* Update trailing submatrix within the band.
333	*
334	IF( JU.GT.J )
335	$ CALL SGER( KM, JU-J, -ONE, AB( KV+2, J ), 1,
336	$ AB( KV, J+1 ), LDAB-1, AB( KV+1, J+1 ),
337	$ LDAB-1 )
338	END IF
339	ELSE
340	*
341	* If pivot is zero, set INFO to the index of the pivot
342	* unless a zero pivot has already been found.
343	*
344	IF( INFO.EQ.0 )
345	$ INFO = J
346	END IF
347	40 CONTINUE
348	RETURN
349	*
350	* End of SGBTF2
351	*
352	END
353
354	C ==================================================================
355
356	SUBROUTINE SGBTRF( M, N, KL, KU, AB, LDAB, IPIV, INFO )
357	*
358	* -- LAPACK routine (version 2.0) --
359	* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
360	* Courant Institute, Argonne National Lab, and Rice University
361	* February 29, 1992
362	*
363	* .. Scalar Arguments ..
364	INTEGER INFO, KL, KU, LDAB, M, N
365	* ..
366	* .. Array Arguments ..
367	INTEGER IPIV( * )
368	REAL AB( LDAB, * )
369	* ..
370	*
371	* Purpose
372	* =======
373	*
374	* SGBTRF computes an LU factorization of a real m-by-n band matrix A
375	* using partial pivoting with row interchanges.
376	*
377	* This is the blocked version of the algorithm, calling Level 3 BLAS.
378	*
379	* Arguments
380	* =========
381	*
382	* M (input) INTEGER
383	* The number of rows of the matrix A. M >= 0.
384	*
385	* N (input) INTEGER
386	* The number of columns of the matrix A. N >= 0.
387	*
388	* KL (input) INTEGER
389	* The number of subdiagonals within the band of A. KL >= 0.
390	*
391	* KU (input) INTEGER
392	* The number of superdiagonals within the band of A. KU >= 0.
393	*
394	* AB (input/output) REAL array, dimension (LDAB,N)
395	* On entry, the matrix A in band storage, in rows KL+1 to
396	* 2*KL+KU+1; rows 1 to KL of the array need not be set.
397	* The j-th column of A is stored in the j-th column of the
398	* array AB as follows:
399	* AB(kl+ku+1+i-j,j) = A(i,j) for max(1,j-ku)<=i<=min(m,j+kl)
400	*
401	* On exit, details of the factorization: U is stored as an
402	* upper triangular band matrix with KL+KU superdiagonals in
403	* rows 1 to KL+KU+1, and the multipliers used during the
404	* factorization are stored in rows KL+KU+2 to 2*KL+KU+1.
405	* See below for further details.
406	*
407	* LDAB (input) INTEGER
408	* The leading dimension of the array AB. LDAB >= 2*KL+KU+1.
409	*
410	* IPIV (output) INTEGER array, dimension (min(M,N))
411	* The pivot indices; for 1 <= i <= min(M,N), row i of the
412	* matrix was interchanged with row IPIV(i).
413	*
414	* INFO (output) INTEGER
415	* = 0: successful exit
416	* < 0: if INFO = -i, the i-th argument had an illegal value
417	* > 0: if INFO = +i, U(i,i) is exactly zero. The factorization
418	* has been completed, but the factor U is exactly
419	* singular, and division by zero will occur if it is used
420	* to solve a system of equations.
421	*
422	* Further Details
423	* ===============
424	*
425	* The band storage scheme is illustrated by the following example, when
426	* M = N = 6, KL = 2, KU = 1:
427	*
428	* On entry: On exit:
429	*
430	* * * * + + + * * * u14 u25 u36
431	* * * + + + + * * u13 u24 u35 u46
432	* * a12 a23 a34 a45 a56 * u12 u23 u34 u45 u56
433	* a11 a22 a33 a44 a55 a66 u11 u22 u33 u44 u55 u66
434	* a21 a32 a43 a54 a65 * m21 m32 m43 m54 m65 *
435	* a31 a42 a53 a64 * * m31 m42 m53 m64 * *
436	*
437	* Array elements marked * are not used by the routine; elements marked
438	* + need not be set on entry, but are required by the routine to store
439	* elements of U because of fill-in resulting from the row interchanges.
440	*
441	* =====================================================================
442	*
443	* .. Parameters ..
444	REAL ONE, ZERO
445	PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 )
446	INTEGER NBMAX, LDWORK
447	PARAMETER ( NBMAX = 64, LDWORK = NBMAX+1 )
448	* ..
449	* .. Local Scalars ..
450	INTEGER I, I2, I3, II, IP, J, J2, J3, JB, JJ, JM, JP,
451	$ JU, K2, KM, KV, NB, NW
452	REAL TEMP
453	* ..
454	* .. Local Arrays ..
455	REAL WORK13( LDWORK, NBMAX ),
456	$ WORK31( LDWORK, NBMAX )
457	* ..
458	* .. External Functions ..
459	INTEGER ILAENV, ISAMAX
460	EXTERNAL ILAENV, ISAMAX
461	* ..
462	* .. External Subroutines ..
463	EXTERNAL SCOPY, SGBTF2, SGEMM, SGER, SLASWP, SSCAL,
464	$ SSWAP, STRSM, XERBLA
465	* ..
466	* .. Intrinsic Functions ..
467	INTRINSIC MAX, MIN
468	* ..
469	* .. Executable Statements ..
470	*
471	* KV is the number of superdiagonals in the factor U, allowing for
472	* fill-in
473	*
474	KV = KU + KL
475	*
476	* Test the input parameters.
477	*
478	INFO = 0
479	IF( M.LT.0 ) THEN
480	INFO = -1
481	ELSE IF( N.LT.0 ) THEN
482	INFO = -2
483	ELSE IF( KL.LT.0 ) THEN
484	INFO = -3
485	ELSE IF( KU.LT.0 ) THEN
486	INFO = -4
487	ELSE IF( LDAB.LT.KL+KV+1 ) THEN
488	INFO = -6
489	END IF
490	IF( INFO.NE.0 ) THEN
491	CALL XERBLA( 'SGBTRF', -INFO )
492	RETURN
493	END IF
494	*
495	* Quick return if possible
496	*
497	IF( M.EQ.0 .OR. N.EQ.0 )
498	$ RETURN
499	*
500	* Determine the block size for this environment
501	*
502	NB = ILAENV( 1, 'SGBTRF', ' ', M, N, KL, KU )
503	*
504	* The block size must not exceed the limit set by the size of the
505	* local arrays WORK13 and WORK31.
506	*
507	NB = MIN( NB, NBMAX )
508	*
509	IF( NB.LE.1 .OR. NB.GT.KL ) THEN
510	*
511	* Use unblocked code
512	*
513	CALL SGBTF2( M, N, KL, KU, AB, LDAB, IPIV, INFO )
514	ELSE
515	*
516	* Use blocked code
517	*
518	* Zero the superdiagonal elements of the work array WORK13
519	*
520	DO 20 J = 1, NB
521	DO 10 I = 1, J - 1
522	WORK13( I, J ) = ZERO
523	10 CONTINUE
524	20 CONTINUE
525	*
526	* Zero the subdiagonal elements of the work array WORK31
527	*
528	DO 40 J = 1, NB
529	DO 30 I = J + 1, NB
530	WORK31( I, J ) = ZERO
531	30 CONTINUE
532	40 CONTINUE
533	*
534	* Gaussian elimination with partial pivoting
535	*
536	* Set fill-in elements in columns KU+2 to KV to zero
537	*
538	DO 60 J = KU + 2, MIN( KV, N )
539	DO 50 I = KV - J + 2, KL
540	AB( I, J ) = ZERO
541	50 CONTINUE
542	60 CONTINUE
543	*
544	* JU is the index of the last column affected by the current
545	* stage of the factorization
546	*
547	JU = 1
548	*
549	DO 180 J = 1, MIN( M, N ), NB
550	JB = MIN( NB, MIN( M, N )-J+1 )
551	*
552	* The active part of the matrix is partitioned
553	*
554	* A11 A12 A13
555	* A21 A22 A23
556	* A31 A32 A33
557	*
558	* Here A11, A21 and A31 denote the current block of JB columns
559	* which is about to be factorized. The number of rows in the
560	* partitioning are JB, I2, I3 respectively, and the numbers
561	* of columns are JB, J2, J3. The superdiagonal elements of A13
562	* and the subdiagonal elements of A31 lie outside the band.
563	*
564	I2 = MIN( KL-JB, M-J-JB+1 )
565	I3 = MIN( JB, M-J-KL+1 )
566	*
567	* J2 and J3 are computed after JU has been updated.
568	*
569	* Factorize the current block of JB columns
570	*
571	DO 80 JJ = J, J + JB - 1
572	*
573	* Set fill-in elements in column JJ+KV to zero
574	*
575	IF( JJ+KV.LE.N ) THEN
576	DO 70 I = 1, KL
577	AB( I, JJ+KV ) = ZERO
578	70 CONTINUE
579	END IF
580	*
581	* Find pivot and test for singularity. KM is the number of
582	* subdiagonal elements in the current column.
583	*
584	KM = MIN( KL, M-JJ )
585	JP = ISAMAX( KM+1, AB( KV+1, JJ ), 1 )
586	IPIV( JJ ) = JP + JJ - J
587	IF( AB( KV+JP, JJ ).NE.ZERO ) THEN
588	JU = MAX( JU, MIN( JJ+KU+JP-1, N ) )
589	IF( JP.NE.1 ) THEN
590	*
591	* Apply interchange to columns J to J+JB-1
592	*
593	IF( JP+JJ-1.LT.J+KL ) THEN
594	*
595	CALL SSWAP( JB, AB( KV+1+JJ-J, J ), LDAB-1,
596	$ AB( KV+JP+JJ-J, J ), LDAB-1 )
597	ELSE
598	*
599	* The interchange affects columns J to JJ-1 of A31
600	* which are stored in the work array WORK31
601	*
602	CALL SSWAP( JJ-J, AB( KV+1+JJ-J, J ), LDAB-1,
603	$ WORK31( JP+JJ-J-KL, 1 ), LDWORK )
604	CALL SSWAP( J+JB-JJ, AB( KV+1, JJ ), LDAB-1,
605	$ AB( KV+JP, JJ ), LDAB-1 )
606	END IF
607	END IF
608	*
609	* Compute multipliers
610	*
611	CALL SSCAL( KM, ONE / AB( KV+1, JJ ), AB( KV+2, JJ ),
612	$ 1 )
613	*
614	* Update trailing submatrix within the band and within
615	* the current block. JM is the index of the last column
616	* which needs to be updated.
617	*
618	JM = MIN( JU, J+JB-1 )
619	IF( JM.GT.JJ )
620	$ CALL SGER( KM, JM-JJ, -ONE, AB( KV+2, JJ ), 1,
621	$ AB( KV, JJ+1 ), LDAB-1,
622	$ AB( KV+1, JJ+1 ), LDAB-1 )
623	ELSE
624	*
625	* If pivot is zero, set INFO to the index of the pivot
626	* unless a zero pivot has already been found.
627	*
628	IF( INFO.EQ.0 )
629	$ INFO = JJ
630	END IF
631	*
632	* Copy current column of A31 into the work array WORK31
633	*
634	NW = MIN( JJ-J+1, I3 )
635	IF( NW.GT.0 )
636	$ CALL SCOPY( NW, AB( KV+KL+1-JJ+J, JJ ), 1,
637	$ WORK31( 1, JJ-J+1 ), 1 )
638	80 CONTINUE
639	IF( J+JB.LE.N ) THEN
640	*
641	* Apply the row interchanges to the other blocks.
642	*
643	J2 = MIN( JU-J+1, KV ) - JB
644	J3 = MAX( 0, JU-J-KV+1 )
645	*
646	* Use SLASWP to apply the row interchanges to A12, A22, and
647	* A32.
648	*
649	CALL SLASWP( J2, AB( KV+1-JB, J+JB ), LDAB-1, 1, JB,
650	$ IPIV( J ), 1 )
651	*
652	* Adjust the pivot indices.
653	*
654	DO 90 I = J, J + JB - 1
655	IPIV( I ) = IPIV( I ) + J - 1
656	90 CONTINUE
657	*
658	* Apply the row interchanges to A13, A23, and A33
659	* columnwise.
660	*
661	K2 = J - 1 + JB + J2
662	DO 110 I = 1, J3
663	JJ = K2 + I
664	DO 100 II = J + I - 1, J + JB - 1
665	IP = IPIV( II )
666	IF( IP.NE.II ) THEN
667	TEMP = AB( KV+1+II-JJ, JJ )
668	AB( KV+1+II-JJ, JJ ) = AB( KV+1+IP-JJ, JJ )
669	AB( KV+1+IP-JJ, JJ ) = TEMP
670	END IF
671	100 CONTINUE
672	110 CONTINUE
673	*
674	* Update the relevant part of the trailing submatrix
675	*
676	IF( J2.GT.0 ) THEN
677	*
678	* Update A12
679	*
680	CALL STRSM( 'Left', 'Lower', 'No transpose', 'Unit',
681	$ JB, J2, ONE, AB( KV+1, J ), LDAB-1,
682	$ AB( KV+1-JB, J+JB ), LDAB-1 )
683	*
684	IF( I2.GT.0 ) THEN
685	*
686	* Update A22
687	*
688	CALL SGEMM( 'No transpose', 'No transpose', I2, J2,
689	$ JB, -ONE, AB( KV+1+JB, J ), LDAB-1,
690	$ AB( KV+1-JB, J+JB ), LDAB-1, ONE,
691	$ AB( KV+1, J+JB ), LDAB-1 )
692	END IF
693	*
694	IF( I3.GT.0 ) THEN
695	*
696	* Update A32
697	*
698	CALL SGEMM( 'No transpose', 'No transpose', I3, J2,
699	$ JB, -ONE, WORK31, LDWORK,
700	$ AB( KV+1-JB, J+JB ), LDAB-1, ONE,
701	$ AB( KV+KL+1-JB, J+JB ), LDAB-1 )
702	END IF
703	END IF
704	*
705	IF( J3.GT.0 ) THEN
706	*
707	* Copy the lower triangle of A13 into the work array
708	* WORK13
709	*
710	DO 130 JJ = 1, J3
711	DO 120 II = JJ, JB
712	WORK13( II, JJ ) = AB( II-JJ+1, JJ+J+KV-1 )
713	120 CONTINUE
714	130 CONTINUE
715	*
716	* Update A13 in the work array
717	*
718	CALL STRSM( 'Left', 'Lower', 'No transpose', 'Unit',
719	$ JB, J3, ONE, AB( KV+1, J ), LDAB-1,
720	$ WORK13, LDWORK )
721	*
722	IF( I2.GT.0 ) THEN
723	*
724	* Update A23
725	*
726	CALL SGEMM( 'No transpose', 'No transpose', I2, J3,
727	$ JB, -ONE, AB( KV+1+JB, J ), LDAB-1,
728	$ WORK13, LDWORK, ONE, AB( 1+JB, J+KV ),
729	$ LDAB-1 )
730	END IF
731	*
732	IF( I3.GT.0 ) THEN
733	*
734	* Update A33
735	*
736	CALL SGEMM( 'No transpose', 'No transpose', I3, J3,
737	$ JB, -ONE, WORK31, LDWORK, WORK13,
738	$ LDWORK, ONE, AB( 1+KL, J+KV ), LDAB-1 )
739	END IF
740	*
741	* Copy the lower triangle of A13 back into place
742	*
743	DO 150 JJ = 1, J3
744	DO 140 II = JJ, JB
745	AB( II-JJ+1, JJ+J+KV-1 ) = WORK13( II, JJ )
746	140 CONTINUE
747	150 CONTINUE
748	END IF
749	ELSE
750	*
751	* Adjust the pivot indices.
752	*
753	DO 160 I = J, J + JB - 1
754	IPIV( I ) = IPIV( I ) + J - 1
755	160 CONTINUE
756	END IF
757	*
758	* Partially undo the interchanges in the current block to
759	* restore the upper triangular form of A31 and copy the upper
760	* triangle of A31 back into place
761	*
762	DO 170 JJ = J + JB - 1, J, -1
763	JP = IPIV( JJ ) - JJ + 1
764	IF( JP.NE.1 ) THEN
765	*
766	* Apply interchange to columns J to JJ-1
767	*
768	IF( JP+JJ-1.LT.J+KL ) THEN
769	*
770	* The interchange does not affect A31
771	*
772	CALL SSWAP( JJ-J, AB( KV+1+JJ-J, J ), LDAB-1,
773	$ AB( KV+JP+JJ-J, J ), LDAB-1 )
774	ELSE
775	*
776	* The interchange does affect A31
777	*
778	CALL SSWAP( JJ-J, AB( KV+1+JJ-J, J ), LDAB-1,
779	$ WORK31( JP+JJ-J-KL, 1 ), LDWORK )
780	END IF
781	END IF
782	*
783	* Copy the current column of A31 back into place
784	*
785	NW = MIN( I3, JJ-J+1 )
786	IF( NW.GT.0 )
787	$ CALL SCOPY( NW, WORK31( 1, JJ-J+1 ), 1,
788	$ AB( KV+KL+1-JJ+J, JJ ), 1 )
789	170 CONTINUE
790	180 CONTINUE
791	END IF
792	*
793	RETURN
794	*
795	* End of SGBTRF
796	*
797	END
798
799	C =======================================================================
800
801	SUBROUTINE SGBTRS( TRANS, N, KL, KU, NRHS, AB, LDAB, IPIV, B, LDB,
802	$ INFO )
803	*
804	* -- LAPACK routine (version 2.0) --
805	* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
806	* Courant Institute, Argonne National Lab, and Rice University
807	* March 31, 1993
808	*
809	* .. Scalar Arguments ..
810	CHARACTER TRANS
811	INTEGER INFO, KL, KU, LDAB, LDB, N, NRHS
812	* ..
813	* .. Array Arguments ..
814	INTEGER IPIV( * )
815	REAL AB( LDAB, * ), B( LDB, * )
816	* ..
817	*
818	* Purpose
819	* =======
820	*
821	* SGBTRS solves a system of linear equations
822	* A * X = B or A' * X = B
823	* with a general band matrix A using the LU factorization computed
824	* by SGBTRF.
825	*
826	* Arguments
827	* =========
828	*
829	* TRANS (input) CHARACTER*1
830	* Specifies the form of the system of equations.
831	* = 'N': A * X = B (No transpose)
832	* = 'T': A'* X = B (Transpose)
833	* = 'C': A'* X = B (Conjugate transpose = Transpose)
834	*
835	* N (input) INTEGER
836	* The order of the matrix A. N >= 0.
837	*
838	* KL (input) INTEGER
839	* The number of subdiagonals within the band of A. KL >= 0.
840	*
841	* KU (input) INTEGER
842	* The number of superdiagonals within the band of A. KU >= 0.
843	*
844	* NRHS (input) INTEGER
845	* The number of right hand sides, i.e., the number of columns
846	* of the matrix B. NRHS >= 0.
847	*
848	* AB (input) REAL array, dimension (LDAB,N)
849	* Details of the LU factorization of the band matrix A, as
850	* computed by SGBTRF. U is stored as an upper triangular band
851	* matrix with KL+KU superdiagonals in rows 1 to KL+KU+1, and
852	* the multipliers used during the factorization are stored in
853	* rows KL+KU+2 to 2*KL+KU+1.
854	*
855	* LDAB (input) INTEGER
856	* The leading dimension of the array AB. LDAB >= 2*KL+KU+1.
857	*
858	* IPIV (input) INTEGER array, dimension (N)
859	* The pivot indices; for 1 <= i <= N, row i of the matrix was
860	* interchanged with row IPIV(i).
861	*
862	* B (input/output) REAL array, dimension (LDB,NRHS)
863	* On entry, the right hand side matrix B.
864	* On exit, the solution matrix X.
865	*
866	* LDB (input) INTEGER
867	* The leading dimension of the array B. LDB >= max(1,N).
868	*
869	* INFO (output) INTEGER
870	* = 0: successful exit
871	* < 0: if INFO = -i, the i-th argument had an illegal value
872	*
873	* =====================================================================
874	*
875	* .. Parameters ..
876	REAL ONE
877	PARAMETER ( ONE = 1.0E+0 )
878	* ..
879	* .. Local Scalars ..
880	LOGICAL LNOTI, NOTRAN
881	INTEGER I, J, KD, L, LM
882	* ..
883	* .. External Functions ..
884	LOGICAL LSAME
885	EXTERNAL LSAME
886	* ..
887	* .. External Subroutines ..
888	EXTERNAL SGEMV, SGER, SSWAP, STBSV, XERBLA
889	* ..
890	* .. Intrinsic Functions ..
891	INTRINSIC MAX, MIN
892	* ..
893	* .. Executable Statements ..
894	*
895	* Test the input parameters.
896	*
897	INFO = 0
898	NOTRAN = LSAME( TRANS, 'N' )
899	IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT.
900	$ LSAME( TRANS, 'C' ) ) THEN
901	INFO = -1
902	ELSE IF( N.LT.0 ) THEN
903	INFO = -2
904	ELSE IF( KL.LT.0 ) THEN
905	INFO = -3
906	ELSE IF( KU.LT.0 ) THEN
907	INFO = -4
908	ELSE IF( NRHS.LT.0 ) THEN
909	INFO = -5
910	ELSE IF( LDAB.LT.( 2*KL+KU+1 ) ) THEN
911	INFO = -7
912	ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
913	INFO = -10
914	END IF
915	IF( INFO.NE.0 ) THEN
916	CALL XERBLA( 'SGBTRS', -INFO )
917	RETURN
918	END IF
919	*
920	* Quick return if possible
921	*
922	IF( N.EQ.0 .OR. NRHS.EQ.0 )
923	$ RETURN
924	*
925	KD = KU + KL + 1
926	LNOTI = KL.GT.0
927	*
928	IF( NOTRAN ) THEN
929	*
930	* Solve A*X = B.
931	*
932	* Solve L*X = B, overwriting B with X.
933	*
934	* L is represented as a product of permutations and unit lower
935	* triangular matrices L = P(1) * L(1) * ... * P(n-1) * L(n-1),
936	* where each transformation L(i) is a rank-one modification of
937	* the identity matrix.
938	*
939	IF( LNOTI ) THEN
940	DO 10 J = 1, N - 1
941	LM = MIN( KL, N-J )
942	L = IPIV( J )
943	IF( L.NE.J )
944	$ CALL SSWAP( NRHS, B( L, 1 ), LDB, B( J, 1 ), LDB )
945	CALL SGER( LM, NRHS, -ONE, AB( KD+1, J ), 1, B( J, 1 ),
946	$ LDB, B( J+1, 1 ), LDB )
947	10 CONTINUE
948	END IF
949	*
950	DO 20 I = 1, NRHS
951	*
952	* Solve U*X = B, overwriting B with X.
953	*
954	CALL STBSV( 'Upper', 'No transpose', 'Non-unit', N, KL+KU,
955	$ AB, LDAB, B( 1, I ), 1 )
956	20 CONTINUE
957	*
958	ELSE
959	*
960	* Solve A'*X = B.
961	*
962	DO 30 I = 1, NRHS
963	*
964	* Solve U'*X = B, overwriting B with X.
965	*
966	CALL STBSV( 'Upper', 'Transpose', 'Non-unit', N, KL+KU, AB,
967	$ LDAB, B( 1, I ), 1 )
968	30 CONTINUE
969	*
970	* Solve L'*X = B, overwriting B with X.
971	*
972	IF( LNOTI ) THEN
973	DO 40 J = N - 1, 1, -1
974	LM = MIN( KL, N-J )
975	CALL SGEMV( 'Transpose', LM, NRHS, -ONE, B( J+1, 1 ),
976	$ LDB, AB( KD+1, J ), 1, ONE, B( J, 1 ), LDB )
977	L = IPIV( J )
978	IF( L.NE.J )
979	$ CALL SSWAP( NRHS, B( L, 1 ), LDB, B( J, 1 ), LDB )
980	40 CONTINUE
981	END IF
982	END IF
983	RETURN
984	*
985	* End of SGBTRS
986	*
987	END
988	C =======================================================================
989	SUBROUTINE SLASWP( N, A, LDA, K1, K2, IPIV, INCX )
990	*
991	* -- LAPACK auxiliary routine (version 2.0) --
992	* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
993	* Courant Institute, Argonne National Lab, and Rice University
994	* October 31, 1992
995	*
996	* .. Scalar Arguments ..
997	INTEGER INCX, K1, K2, LDA, N
998	* ..
999	* .. Array Arguments ..
1000	INTEGER IPIV( * )
1001	REAL A( LDA, * )
1002	* ..
1003	*
1004	* Purpose
1005	* =======
1006	*
1007	* SLASWP performs a series of row interchanges on the matrix A.
1008	* One row interchange is initiated for each of rows K1 through K2 of A.
1009	*
1010	* Arguments
1011	* =========
1012	*
1013	* N (input) INTEGER
1014	* The number of columns of the matrix A.
1015	*
1016	* A (input/output) REAL array, dimension (LDA,N)
1017	* On entry, the matrix of column dimension N to which the row
1018	* interchanges will be applied.
1019	* On exit, the permuted matrix.
1020	*
1021	* LDA (input) INTEGER
1022	* The leading dimension of the array A.
1023	*
1024	* K1 (input) INTEGER
1025	* The first element of IPIV for which a row interchange will
1026	* be done.
1027	*
1028	* K2 (input) INTEGER
1029	* The last element of IPIV for which a row interchange will
1030	* be done.
1031	*
1032	* IPIV (input) INTEGER array, dimension (M*abs(INCX))
1033	* The vector of pivot indices. Only the elements in positions
1034	* K1 through K2 of IPIV are accessed.
1035	* IPIV(K) = L implies rows K and L are to be interchanged.
1036	*
1037	* INCX (input) INTEGER
1038	* The increment between successive values of IPIV. If IPIV
1039	* is negative, the pivots are applied in reverse order.
1040	*
1041	* =====================================================================
1042	*
1043	* .. Local Scalars ..
1044	INTEGER I, IP, IX
1045	* ..
1046	* .. External Subroutines ..
1047	EXTERNAL SSWAP
1048	* ..
1049	* .. Executable Statements ..
1050	*
1051	* Interchange row I with row IPIV(I) for each of rows K1 through K2.
1052	*
1053	IF( INCX.EQ.0 )
1054	$ RETURN
1055	IF( INCX.GT.0 ) THEN
1056	IX = K1
1057	ELSE
1058	IX = 1 + ( 1-K2 )*INCX
1059	END IF
1060	IF( INCX.EQ.1 ) THEN
1061	DO 10 I = K1, K2
1062	IP = IPIV( I )
1063	IF( IP.NE.I )
1064	$ CALL SSWAP( N, A( I, 1 ), LDA, A( IP, 1 ), LDA )
1065	10 CONTINUE
1066	ELSE IF( INCX.GT.1 ) THEN
1067	DO 20 I = K1, K2
1068	IP = IPIV( IX )
1069	IF( IP.NE.I )
1070	$ CALL SSWAP( N, A( I, 1 ), LDA, A( IP, 1 ), LDA )
1071	IX = IX + INCX
1072	20 CONTINUE
1073	ELSE IF( INCX.LT.0 ) THEN
1074	DO 30 I = K2, K1, -1
1075	IP = IPIV( IX )
1076	IF( IP.NE.I )
1077	$ CALL SSWAP( N, A( I, 1 ), LDA, A( IP, 1 ), LDA )
1078	IX = IX + INCX
1079	30 CONTINUE
1080	END IF
1081	*
1082	RETURN
1083	*
1084	* End of SLASWP
1085	*
1086	END
1087	C =====================================================================
1088	INTEGER FUNCTION ILAENV( ISPEC, NAME, OPTS, N1, N2, N3,
1089	$ N4 )
1090	*
1091	* -- LAPACK auxiliary routine (version 2.0) --
1092	* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
1093	* Courant Institute, Argonne National Lab, and Rice University
1094	* September 30, 1994
1095	*
1096	* .. Scalar Arguments ..
1097	CHARACTER( ) NAME, OPTS
1098	INTEGER ISPEC, N1, N2, N3, N4
1099	* ..
1100	*
1101	* Purpose
1102	* =======
1103	*
1104	* ILAENV is called from the LAPACK routines to choose problem-dependent
1105	* parameters for the local environment. See ISPEC for a description of
1106	* the parameters.
1107	*
1108	* This version provides a set of parameters which should give good,
1109	* but not optimal, performance on many of the currently available
1110	* computers. Users are encouraged to modify this subroutine to set
1111	* the tuning parameters for their particular machine using the option
1112	* and problem size information in the arguments.
1113	*
1114	* This routine will not function correctly if it is converted to all
1115	* lower case. Converting it to all upper case is allowed.
1116	*
1117	* Arguments
1118	* =========
1119	*
1120	* ISPEC (input) INTEGER
1121	* Specifies the parameter to be returned as the value of
1122	* ILAENV.
1123	* = 1: the optimal blocksize; if this value is 1, an unblocked
1124	* algorithm will give the best performance.
1125	* = 2: the minimum block size for which the block routine
1126	* should be used; if the usable block size is less than
1127	* this value, an unblocked routine should be used.
1128	* = 3: the crossover point (in a block routine, for N less
1129	* than this value, an unblocked routine should be used)
1130	* = 4: the number of shifts, used in the nonsymmetric
1131	* eigenvalue routines
1132	* = 5: the minimum column dimension for blocking to be used;
1133	* rectangular blocks must have dimension at least k by m,
1134	* where k is given by ILAENV(2,...) and m by ILAENV(5,...)
1135	* = 6: the crossover point for the SVD (when reducing an m by n
1136	* matrix to bidiagonal form, if max(m,n)/min(m,n) exceeds
1137	* this value, a QR factorization is used first to reduce
1138	* the matrix to a triangular form.)
1139	* = 7: the number of processors
1140	* = 8: the crossover point for the multishift QR and QZ methods
1141	* for nonsymmetric eigenvalue problems.
1142	*
1143	* NAME (input) CHARACTER()
1144	* The name of the calling subroutine, in either upper case or
1145	* lower case.
1146	*
1147	* OPTS (input) CHARACTER()
1148	* The character options to the subroutine NAME, concatenated
1149	* into a single character string. For example, UPLO = 'U',
1150	* TRANS = 'T', and DIAG = 'N' for a triangular routine would
1151	* be specified as OPTS = 'UTN'.
1152	*
1153	* N1 (input) INTEGER
1154	* N2 (input) INTEGER
1155	* N3 (input) INTEGER
1156	* N4 (input) INTEGER
1157	* Problem dimensions for the subroutine NAME; these may not all
1158	* be required.
1159	*
1160	* (ILAENV) (output) INTEGER
1161	* >= 0: the value of the parameter specified by ISPEC
1162	* < 0: if ILAENV = -k, the k-th argument had an illegal value.
1163	*
1164	* Further Details
1165	* ===============
1166	*
1167	* The following conventions have been used when calling ILAENV from the
1168	* LAPACK routines:
1169	* 1) OPTS is a concatenation of all of the character options to
1170	* subroutine NAME, in the same order that they appear in the
1171	* argument list for NAME, even if they are not used in determining
1172	* the value of the parameter specified by ISPEC.
1173	* 2) The problem dimensions N1, N2, N3, N4 are specified in the order
1174	* that they appear in the argument list for NAME. N1 is used
1175	* first, N2 second, and so on, and unused problem dimensions are
1176	* passed a value of -1.
1177	* 3) The parameter value returned by ILAENV is checked for validity in
1178	* the calling subroutine. For example, ILAENV is used to retrieve
1179	* the optimal blocksize for STRTRI as follows:
1180	*
1181	* NB = ILAENV( 1, 'STRTRI', UPLO // DIAG, N, -1, -1, -1 )
1182	* IF( NB.LE.1 ) NB = MAX( 1, N )
1183	*
1184	* =====================================================================
1185	*
1186	* .. Local Scalars ..
1187	LOGICAL CNAME, SNAME
1188	CHARACTER*1 C1
1189	CHARACTER*2 C2, C4
1190	CHARACTER*3 C3
1191	CHARACTER*6 SUBNAM
1192	INTEGER I, IC, IZ, NB, NBMIN, NX
1193	* ..
1194	* .. Intrinsic Functions ..
1195	INTRINSIC CHAR, ICHAR, INT, MIN, REAL
1196	* ..
1197	* .. Executable Statements ..
1198	*
1199	GO TO ( 100, 100, 100, 400, 500, 600, 700, 800 ) ISPEC
1200	*
1201	* Invalid value for ISPEC
1202	*
1203	ILAENV = -1
1204	RETURN
1205	*
1206	100 CONTINUE
1207	*
1208	* Convert NAME to upper case if the first character is lower case.
1209	*
1210	ILAENV = 1
1211	SUBNAM = NAME
1212	IC = ICHAR( SUBNAM( 1:1 ) )
1213	IZ = ICHAR( 'Z' )
1214	IF( IZ.EQ.90 .OR. IZ.EQ.122 ) THEN
1215	*
1216	* ASCII character set
1217	*
1218	IF( IC.GE.97 .AND. IC.LE.122 ) THEN
1219	SUBNAM( 1:1 ) = CHAR( IC-32 )
1220	DO 10 I = 2, 6
1221	IC = ICHAR( SUBNAM( I:I ) )
1222	IF( IC.GE.97 .AND. IC.LE.122 )
1223	$ SUBNAM( I:I ) = CHAR( IC-32 )
1224	10 CONTINUE
1225	END IF
1226	*
1227	ELSE IF( IZ.EQ.233 .OR. IZ.EQ.169 ) THEN
1228	*
1229	* EBCDIC character set
1230	*
1231	IF( ( IC.GE.129 .AND. IC.LE.137 ) .OR.
1232	$ ( IC.GE.145 .AND. IC.LE.153 ) .OR.
1233	$ ( IC.GE.162 .AND. IC.LE.169 ) ) THEN
1234	SUBNAM( 1:1 ) = CHAR( IC+64 )
1235	DO 20 I = 2, 6
1236	IC = ICHAR( SUBNAM( I:I ) )
1237	IF( ( IC.GE.129 .AND. IC.LE.137 ) .OR.
1238	$ ( IC.GE.145 .AND. IC.LE.153 ) .OR.
1239	$ ( IC.GE.162 .AND. IC.LE.169 ) )
1240	$ SUBNAM( I:I ) = CHAR( IC+64 )
1241	20 CONTINUE
1242	END IF
1243	*
1244	ELSE IF( IZ.EQ.218 .OR. IZ.EQ.250 ) THEN
1245	*
1246	* Prime machines: ASCII+128
1247	*
1248	IF( IC.GE.225 .AND. IC.LE.250 ) THEN
1249	SUBNAM( 1:1 ) = CHAR( IC-32 )
1250	DO 30 I = 2, 6
1251	IC = ICHAR( SUBNAM( I:I ) )
1252	IF( IC.GE.225 .AND. IC.LE.250 )
1253	$ SUBNAM( I:I ) = CHAR( IC-32 )
1254	30 CONTINUE
1255	END IF
1256	END IF
1257	*
1258	C1 = SUBNAM( 1:1 )
1259	SNAME = C1.EQ.'S' .OR. C1.EQ.'D'
1260	CNAME = C1.EQ.'C' .OR. C1.EQ.'Z'
1261	IF( .NOT.( CNAME .OR. SNAME ) )
1262	$ RETURN
1263	C2 = SUBNAM( 2:3 )
1264	C3 = SUBNAM( 4:6 )
1265	C4 = C3( 2:3 )
1266	*
1267	GO TO ( 110, 200, 300 ) ISPEC
1268	*
1269	110 CONTINUE
1270	*
1271	* ISPEC = 1: block size
1272	*
1273	* In these examples, separate code is provided for setting NB for
1274	* real and complex. We assume that NB will take the same value in
1275	* single or double precision.
1276	*
1277	NB = 1
1278	*
1279	IF( C2.EQ.'GE' ) THEN
1280	IF( C3.EQ.'TRF' ) THEN
1281	IF( SNAME ) THEN
1282	NB = 64
1283	ELSE
1284	NB = 64
1285	END IF
1286	ELSE IF( C3.EQ.'QRF' .OR. C3.EQ.'RQF' .OR. C3.EQ.'LQF' .OR.
1287	$ C3.EQ.'QLF' ) THEN
1288	IF( SNAME ) THEN
1289	NB = 32
1290	ELSE
1291	NB = 32
1292	END IF
1293	ELSE IF( C3.EQ.'HRD' ) THEN
1294	IF( SNAME ) THEN
1295	NB = 32
1296	ELSE
1297	NB = 32
1298	END IF
1299	ELSE IF( C3.EQ.'BRD' ) THEN
1300	IF( SNAME ) THEN
1301	NB = 32
1302	ELSE
1303	NB = 32
1304	END IF
1305	ELSE IF( C3.EQ.'TRI' ) THEN
1306	IF( SNAME ) THEN
1307	NB = 64
1308	ELSE
1309	NB = 64
1310	END IF
1311	END IF
1312	ELSE IF( C2.EQ.'PO' ) THEN
1313	IF( C3.EQ.'TRF' ) THEN
1314	IF( SNAME ) THEN
1315	NB = 64
1316	ELSE
1317	NB = 64
1318	END IF
1319	END IF
1320	ELSE IF( C2.EQ.'SY' ) THEN
1321	IF( C3.EQ.'TRF' ) THEN
1322	IF( SNAME ) THEN
1323	NB = 64
1324	ELSE
1325	NB = 64
1326	END IF
1327	ELSE IF( SNAME .AND. C3.EQ.'TRD' ) THEN
1328	NB = 1
1329	ELSE IF( SNAME .AND. C3.EQ.'GST' ) THEN
1330	NB = 64
1331	END IF
1332	ELSE IF( CNAME .AND. C2.EQ.'HE' ) THEN
1333	IF( C3.EQ.'TRF' ) THEN
1334	NB = 64
1335	ELSE IF( C3.EQ.'TRD' ) THEN
1336	NB = 1
1337	ELSE IF( C3.EQ.'GST' ) THEN
1338	NB = 64
1339	END IF
1340	ELSE IF( SNAME .AND. C2.EQ.'OR' ) THEN
1341	IF( C3( 1:1 ).EQ.'G' ) THEN
1342	IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR.
1343	$ C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR.
1344	$ C4.EQ.'BR' ) THEN
1345	NB = 32
1346	END IF
1347	ELSE IF( C3( 1:1 ).EQ.'M' ) THEN
1348	IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR.
1349	$ C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR.
1350	$ C4.EQ.'BR' ) THEN
1351	NB = 32
1352	END IF
1353	END IF
1354	ELSE IF( CNAME .AND. C2.EQ.'UN' ) THEN
1355	IF( C3( 1:1 ).EQ.'G' ) THEN
1356	IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR.
1357	$ C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR.
1358	$ C4.EQ.'BR' ) THEN
1359	NB = 32
1360	END IF
1361	ELSE IF( C3( 1:1 ).EQ.'M' ) THEN
1362	IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR.
1363	$ C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR.
1364	$ C4.EQ.'BR' ) THEN
1365	NB = 32
1366	END IF
1367	END IF
1368	ELSE IF( C2.EQ.'GB' ) THEN
1369	IF( C3.EQ.'TRF' ) THEN
1370	IF( SNAME ) THEN
1371	IF( N4.LE.64 ) THEN
1372	NB = 1
1373	ELSE
1374	NB = 32
1375	END IF
1376	ELSE
1377	IF( N4.LE.64 ) THEN
1378	NB = 1
1379	ELSE
1380	NB = 32
1381	END IF
1382	END IF
1383	END IF
1384	ELSE IF( C2.EQ.'PB' ) THEN
1385	IF( C3.EQ.'TRF' ) THEN
1386	IF( SNAME ) THEN
1387	IF( N2.LE.64 ) THEN
1388	NB = 1
1389	ELSE
1390	NB = 32
1391	END IF
1392	ELSE
1393	IF( N2.LE.64 ) THEN
1394	NB = 1
1395	ELSE
1396	NB = 32
1397	END IF
1398	END IF
1399	END IF
1400	ELSE IF( C2.EQ.'TR' ) THEN
1401	IF( C3.EQ.'TRI' ) THEN
1402	IF( SNAME ) THEN
1403	NB = 64
1404	ELSE
1405	NB = 64
1406	END IF
1407	END IF
1408	ELSE IF( C2.EQ.'LA' ) THEN
1409	IF( C3.EQ.'UUM' ) THEN
1410	IF( SNAME ) THEN
1411	NB = 64
1412	ELSE
1413	NB = 64
1414	END IF
1415	END IF
1416	ELSE IF( SNAME .AND. C2.EQ.'ST' ) THEN
1417	IF( C3.EQ.'EBZ' ) THEN
1418	NB = 1
1419	END IF
1420	END IF
1421	ILAENV = NB
1422	RETURN
1423	*
1424	200 CONTINUE
1425	*
1426	* ISPEC = 2: minimum block size
1427	*
1428	NBMIN = 2
1429	IF( C2.EQ.'GE' ) THEN
1430	IF( C3.EQ.'QRF' .OR. C3.EQ.'RQF' .OR. C3.EQ.'LQF' .OR.
1431	$ C3.EQ.'QLF' ) THEN
1432	IF( SNAME ) THEN
1433	NBMIN = 2
1434	ELSE
1435	NBMIN = 2
1436	END IF
1437	ELSE IF( C3.EQ.'HRD' ) THEN
1438	IF( SNAME ) THEN
1439	NBMIN = 2
1440	ELSE
1441	NBMIN = 2
1442	END IF
1443	ELSE IF( C3.EQ.'BRD' ) THEN
1444	IF( SNAME ) THEN
1445	NBMIN = 2
1446	ELSE
1447	NBMIN = 2
1448	END IF
1449	ELSE IF( C3.EQ.'TRI' ) THEN
1450	IF( SNAME ) THEN
1451	NBMIN = 2
1452	ELSE
1453	NBMIN = 2
1454	END IF
1455	END IF
1456	ELSE IF( C2.EQ.'SY' ) THEN
1457	IF( C3.EQ.'TRF' ) THEN
1458	IF( SNAME ) THEN
1459	NBMIN = 8
1460	ELSE
1461	NBMIN = 8
1462	END IF
1463	ELSE IF( SNAME .AND. C3.EQ.'TRD' ) THEN
1464	NBMIN = 2
1465	END IF
1466	ELSE IF( CNAME .AND. C2.EQ.'HE' ) THEN
1467	IF( C3.EQ.'TRD' ) THEN
1468	NBMIN = 2
1469	END IF
1470	ELSE IF( SNAME .AND. C2.EQ.'OR' ) THEN
1471	IF( C3( 1:1 ).EQ.'G' ) THEN
1472	IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR.
1473	$ C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR.
1474	$ C4.EQ.'BR' ) THEN
1475	NBMIN = 2
1476	END IF
1477	ELSE IF( C3( 1:1 ).EQ.'M' ) THEN
1478	IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR.
1479	$ C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR.
1480	$ C4.EQ.'BR' ) THEN
1481	NBMIN = 2
1482	END IF
1483	END IF
1484	ELSE IF( CNAME .AND. C2.EQ.'UN' ) THEN
1485	IF( C3( 1:1 ).EQ.'G' ) THEN
1486	IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR.
1487	$ C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR.
1488	$ C4.EQ.'BR' ) THEN
1489	NBMIN = 2
1490	END IF
1491	ELSE IF( C3( 1:1 ).EQ.'M' ) THEN
1492	IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR.
1493	$ C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR.
1494	$ C4.EQ.'BR' ) THEN
1495	NBMIN = 2
1496	END IF
1497	END IF
1498	END IF
1499	ILAENV = NBMIN
1500	RETURN
1501	*
1502	300 CONTINUE
1503	*
1504	* ISPEC = 3: crossover point
1505	*
1506	NX = 0
1507	IF( C2.EQ.'GE' ) THEN
1508	IF( C3.EQ.'QRF' .OR. C3.EQ.'RQF' .OR. C3.EQ.'LQF' .OR.
1509	$ C3.EQ.'QLF' ) THEN
1510	IF( SNAME ) THEN
1511	NX = 128
1512	ELSE
1513	NX = 128
1514	END IF
1515	ELSE IF( C3.EQ.'HRD' ) THEN
1516	IF( SNAME ) THEN
1517	NX = 128
1518	ELSE
1519	NX = 128
1520	END IF
1521	ELSE IF( C3.EQ.'BRD' ) THEN
1522	IF( SNAME ) THEN
1523	NX = 128
1524	ELSE
1525	NX = 128
1526	END IF
1527	END IF
1528	ELSE IF( C2.EQ.'SY' ) THEN
1529	IF( SNAME .AND. C3.EQ.'TRD' ) THEN
1530	NX = 1
1531	END IF
1532	ELSE IF( CNAME .AND. C2.EQ.'HE' ) THEN
1533	IF( C3.EQ.'TRD' ) THEN
1534	NX = 1
1535	END IF
1536	ELSE IF( SNAME .AND. C2.EQ.'OR' ) THEN
1537	IF( C3( 1:1 ).EQ.'G' ) THEN
1538	IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR.
1539	$ C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR.
1540	$ C4.EQ.'BR' ) THEN
1541	NX = 128
1542	END IF
1543	END IF
1544	ELSE IF( CNAME .AND. C2.EQ.'UN' ) THEN
1545	IF( C3( 1:1 ).EQ.'G' ) THEN
1546	IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR.
1547	$ C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR.
1548	$ C4.EQ.'BR' ) THEN
1549	NX = 128
1550	END IF
1551	END IF
1552	END IF
1553	ILAENV = NX
1554	RETURN
1555	*
1556	400 CONTINUE
1557	*
1558	* ISPEC = 4: number of shifts (used by xHSEQR)
1559	*
1560	ILAENV = 6
1561	RETURN
1562	*
1563	500 CONTINUE
1564	*
1565	* ISPEC = 5: minimum column dimension (not used)
1566	*
1567	ILAENV = 2
1568	RETURN
1569	*
1570	600 CONTINUE
1571	*
1572	* ISPEC = 6: crossover point for SVD (used by xGELSS and xGESVD)
1573	*
1574	ILAENV = INT( REAL( MIN( N1, N2 ) )*1.6E0 )
1575	RETURN
1576	*
1577	700 CONTINUE
1578	*
1579	* ISPEC = 7: number of processors (not used)
1580	*
1581	ILAENV = 1
1582	RETURN
1583	*
1584	800 CONTINUE
1585	*
1586	* ISPEC = 8: crossover point for multishift (used by xHSEQR)
1587	*
1588	ILAENV = 50
1589	RETURN
1590	*
1591	* End of ILAENV
1592	*
1593	END
1594

Note: See TracBrowser for help on using the repository browser.

Download in other formats: