1: *> \brief <b> ZGGEV3 computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices (blocked algorithm)</b>
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download ZGGEV3 + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zggev3.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zggev3.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zggev3.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE ZGGEV3( JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHA, BETA,
22: * VL, LDVL, VR, LDVR, WORK, LWORK, RWORK, INFO )
23: *
24: * .. Scalar Arguments ..
25: * CHARACTER JOBVL, JOBVR
26: * INTEGER INFO, LDA, LDB, LDVL, LDVR, LWORK, N
27: * ..
28: * .. Array Arguments ..
29: * DOUBLE PRECISION RWORK( * )
30: * COMPLEX*16 A( LDA, * ), ALPHA( * ), B( LDB, * ),
31: * $ BETA( * ), VL( LDVL, * ), VR( LDVR, * ),
32: * $ WORK( * )
33: * ..
34: *
35: *
36: *> \par Purpose:
37: * =============
38: *>
39: *> \verbatim
40: *>
41: *> ZGGEV3 computes for a pair of N-by-N complex nonsymmetric matrices
42: *> (A,B), the generalized eigenvalues, and optionally, the left and/or
43: *> right generalized eigenvectors.
44: *>
45: *> A generalized eigenvalue for a pair of matrices (A,B) is a scalar
46: *> lambda or a ratio alpha/beta = lambda, such that A - lambda*B is
47: *> singular. It is usually represented as the pair (alpha,beta), as
48: *> there is a reasonable interpretation for beta=0, and even for both
49: *> being zero.
50: *>
51: *> The right generalized eigenvector v(j) corresponding to the
52: *> generalized eigenvalue lambda(j) of (A,B) satisfies
53: *>
54: *> A * v(j) = lambda(j) * B * v(j).
55: *>
56: *> The left generalized eigenvector u(j) corresponding to the
57: *> generalized eigenvalues lambda(j) of (A,B) satisfies
58: *>
59: *> u(j)**H * A = lambda(j) * u(j)**H * B
60: *>
61: *> where u(j)**H is the conjugate-transpose of u(j).
62: *> \endverbatim
63: *
64: * Arguments:
65: * ==========
66: *
67: *> \param[in] JOBVL
68: *> \verbatim
69: *> JOBVL is CHARACTER*1
70: *> = 'N': do not compute the left generalized eigenvectors;
71: *> = 'V': compute the left generalized eigenvectors.
72: *> \endverbatim
73: *>
74: *> \param[in] JOBVR
75: *> \verbatim
76: *> JOBVR is CHARACTER*1
77: *> = 'N': do not compute the right generalized eigenvectors;
78: *> = 'V': compute the right generalized eigenvectors.
79: *> \endverbatim
80: *>
81: *> \param[in] N
82: *> \verbatim
83: *> N is INTEGER
84: *> The order of the matrices A, B, VL, and VR. N >= 0.
85: *> \endverbatim
86: *>
87: *> \param[in,out] A
88: *> \verbatim
89: *> A is COMPLEX*16 array, dimension (LDA, N)
90: *> On entry, the matrix A in the pair (A,B).
91: *> On exit, A has been overwritten.
92: *> \endverbatim
93: *>
94: *> \param[in] LDA
95: *> \verbatim
96: *> LDA is INTEGER
97: *> The leading dimension of A. LDA >= max(1,N).
98: *> \endverbatim
99: *>
100: *> \param[in,out] B
101: *> \verbatim
102: *> B is COMPLEX*16 array, dimension (LDB, N)
103: *> On entry, the matrix B in the pair (A,B).
104: *> On exit, B has been overwritten.
105: *> \endverbatim
106: *>
107: *> \param[in] LDB
108: *> \verbatim
109: *> LDB is INTEGER
110: *> The leading dimension of B. LDB >= max(1,N).
111: *> \endverbatim
112: *>
113: *> \param[out] ALPHA
114: *> \verbatim
115: *> ALPHA is COMPLEX*16 array, dimension (N)
116: *> \endverbatim
117: *>
118: *> \param[out] BETA
119: *> \verbatim
120: *> BETA is COMPLEX*16 array, dimension (N)
121: *> On exit, ALPHA(j)/BETA(j), j=1,...,N, will be the
122: *> generalized eigenvalues.
123: *>
124: *> Note: the quotients ALPHA(j)/BETA(j) may easily over- or
125: *> underflow, and BETA(j) may even be zero. Thus, the user
126: *> should avoid naively computing the ratio alpha/beta.
127: *> However, ALPHA will be always less than and usually
128: *> comparable with norm(A) in magnitude, and BETA always less
129: *> than and usually comparable with norm(B).
130: *> \endverbatim
131: *>
132: *> \param[out] VL
133: *> \verbatim
134: *> VL is COMPLEX*16 array, dimension (LDVL,N)
135: *> If JOBVL = 'V', the left generalized eigenvectors u(j) are
136: *> stored one after another in the columns of VL, in the same
137: *> order as their eigenvalues.
138: *> Each eigenvector is scaled so the largest component has
139: *> abs(real part) + abs(imag. part) = 1.
140: *> Not referenced if JOBVL = 'N'.
141: *> \endverbatim
142: *>
143: *> \param[in] LDVL
144: *> \verbatim
145: *> LDVL is INTEGER
146: *> The leading dimension of the matrix VL. LDVL >= 1, and
147: *> if JOBVL = 'V', LDVL >= N.
148: *> \endverbatim
149: *>
150: *> \param[out] VR
151: *> \verbatim
152: *> VR is COMPLEX*16 array, dimension (LDVR,N)
153: *> If JOBVR = 'V', the right generalized eigenvectors v(j) are
154: *> stored one after another in the columns of VR, in the same
155: *> order as their eigenvalues.
156: *> Each eigenvector is scaled so the largest component has
157: *> abs(real part) + abs(imag. part) = 1.
158: *> Not referenced if JOBVR = 'N'.
159: *> \endverbatim
160: *>
161: *> \param[in] LDVR
162: *> \verbatim
163: *> LDVR is INTEGER
164: *> The leading dimension of the matrix VR. LDVR >= 1, and
165: *> if JOBVR = 'V', LDVR >= N.
166: *> \endverbatim
167: *>
168: *> \param[out] WORK
169: *> \verbatim
170: *> WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
171: *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
172: *> \endverbatim
173: *>
174: *> \param[in] LWORK
175: *> \verbatim
176: *> LWORK is INTEGER
177: *> The dimension of the array WORK.
178: *>
179: *> If LWORK = -1, then a workspace query is assumed; the routine
180: *> only calculates the optimal size of the WORK array, returns
181: *> this value as the first entry of the WORK array, and no error
182: *> message related to LWORK is issued by XERBLA.
183: *> \endverbatim
184: *>
185: *> \param[out] RWORK
186: *> \verbatim
187: *> RWORK is DOUBLE PRECISION array, dimension (8*N)
188: *> \endverbatim
189: *>
190: *> \param[out] INFO
191: *> \verbatim
192: *> INFO is INTEGER
193: *> = 0: successful exit
194: *> < 0: if INFO = -i, the i-th argument had an illegal value.
195: *> =1,...,N:
196: *> The QZ iteration failed. No eigenvectors have been
197: *> calculated, but ALPHA(j) and BETA(j) should be
198: *> correct for j=INFO+1,...,N.
199: *> > N: =N+1: other then QZ iteration failed in DHGEQZ,
200: *> =N+2: error return from DTGEVC.
201: *> \endverbatim
202: *
203: * Authors:
204: * ========
205: *
206: *> \author Univ. of Tennessee
207: *> \author Univ. of California Berkeley
208: *> \author Univ. of Colorado Denver
209: *> \author NAG Ltd.
210: *
211: *> \date January 2015
212: *
213: *> \ingroup complex16GEeigen
214: *
215: * =====================================================================
216: SUBROUTINE ZGGEV3( JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHA, BETA,
217: $ VL, LDVL, VR, LDVR, WORK, LWORK, RWORK, INFO )
218: *
219: * -- LAPACK driver routine (version 3.6.1) --
220: * -- LAPACK is a software package provided by Univ. of Tennessee, --
221: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
222: * January 2015
223: *
224: * .. Scalar Arguments ..
225: CHARACTER JOBVL, JOBVR
226: INTEGER INFO, LDA, LDB, LDVL, LDVR, LWORK, N
227: * ..
228: * .. Array Arguments ..
229: DOUBLE PRECISION RWORK( * )
230: COMPLEX*16 A( LDA, * ), ALPHA( * ), B( LDB, * ),
231: $ BETA( * ), VL( LDVL, * ), VR( LDVR, * ),
232: $ WORK( * )
233: * ..
234: *
235: * =====================================================================
236: *
237: * .. Parameters ..
238: DOUBLE PRECISION ZERO, ONE
239: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
240: COMPLEX*16 CZERO, CONE
241: PARAMETER ( CZERO = ( 0.0D0, 0.0D0 ),
242: $ CONE = ( 1.0D0, 0.0D0 ) )
243: * ..
244: * .. Local Scalars ..
245: LOGICAL ILASCL, ILBSCL, ILV, ILVL, ILVR, LQUERY
246: CHARACTER CHTEMP
247: INTEGER ICOLS, IERR, IHI, IJOBVL, IJOBVR, ILEFT, ILO,
248: $ IN, IRIGHT, IROWS, IRWRK, ITAU, IWRK, JC, JR,
249: $ LWKOPT
250: DOUBLE PRECISION ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS,
251: $ SMLNUM, TEMP
252: COMPLEX*16 X
253: * ..
254: * .. Local Arrays ..
255: LOGICAL LDUMMA( 1 )
256: * ..
257: * .. External Subroutines ..
258: EXTERNAL DLABAD, XERBLA, ZGEQRF, ZGGBAK, ZGGBAL, ZGGHD3,
259: $ ZHGEQZ, ZLACPY, ZLASCL, ZLASET, ZTGEVC, ZUNGQR,
260: $ ZUNMQR
261: * ..
262: * .. External Functions ..
263: LOGICAL LSAME
264: DOUBLE PRECISION DLAMCH, ZLANGE
265: EXTERNAL LSAME, DLAMCH, ZLANGE
266: * ..
267: * .. Intrinsic Functions ..
268: INTRINSIC ABS, DBLE, DIMAG, MAX, SQRT
269: * ..
270: * .. Statement Functions ..
271: DOUBLE PRECISION ABS1
272: * ..
273: * .. Statement Function definitions ..
274: ABS1( X ) = ABS( DBLE( X ) ) + ABS( DIMAG( X ) )
275: * ..
276: * .. Executable Statements ..
277: *
278: * Decode the input arguments
279: *
280: IF( LSAME( JOBVL, 'N' ) ) THEN
281: IJOBVL = 1
282: ILVL = .FALSE.
283: ELSE IF( LSAME( JOBVL, 'V' ) ) THEN
284: IJOBVL = 2
285: ILVL = .TRUE.
286: ELSE
287: IJOBVL = -1
288: ILVL = .FALSE.
289: END IF
290: *
291: IF( LSAME( JOBVR, 'N' ) ) THEN
292: IJOBVR = 1
293: ILVR = .FALSE.
294: ELSE IF( LSAME( JOBVR, 'V' ) ) THEN
295: IJOBVR = 2
296: ILVR = .TRUE.
297: ELSE
298: IJOBVR = -1
299: ILVR = .FALSE.
300: END IF
301: ILV = ILVL .OR. ILVR
302: *
303: * Test the input arguments
304: *
305: INFO = 0
306: LQUERY = ( LWORK.EQ.-1 )
307: IF( IJOBVL.LE.0 ) THEN
308: INFO = -1
309: ELSE IF( IJOBVR.LE.0 ) THEN
310: INFO = -2
311: ELSE IF( N.LT.0 ) THEN
312: INFO = -3
313: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
314: INFO = -5
315: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
316: INFO = -7
317: ELSE IF( LDVL.LT.1 .OR. ( ILVL .AND. LDVL.LT.N ) ) THEN
318: INFO = -11
319: ELSE IF( LDVR.LT.1 .OR. ( ILVR .AND. LDVR.LT.N ) ) THEN
320: INFO = -13
321: ELSE IF( LWORK.LT.MAX( 1, 2*N ) .AND. .NOT.LQUERY ) THEN
322: INFO = -15
323: END IF
324: *
325: * Compute workspace
326: *
327: IF( INFO.EQ.0 ) THEN
328: CALL ZGEQRF( N, N, B, LDB, WORK, WORK, -1, IERR )
329: LWKOPT = MAX( 1, N+INT( WORK( 1 ) ) )
330: CALL ZUNMQR( 'L', 'C', N, N, N, B, LDB, WORK, A, LDA, WORK,
331: $ -1, IERR )
332: LWKOPT = MAX( LWKOPT, N+INT( WORK( 1 ) ) )
333: IF( ILVL ) THEN
334: CALL ZUNGQR( N, N, N, VL, LDVL, WORK, WORK, -1, IERR )
335: LWKOPT = MAX( LWKOPT, N+INT( WORK( 1 ) ) )
336: END IF
337: IF( ILV ) THEN
338: CALL ZGGHD3( JOBVL, JOBVR, N, 1, N, A, LDA, B, LDB, VL,
339: $ LDVL, VR, LDVR, WORK, -1, IERR )
340: LWKOPT = MAX( LWKOPT, N+INT( WORK( 1 ) ) )
341: CALL ZHGEQZ( 'S', JOBVL, JOBVR, N, 1, N, A, LDA, B, LDB,
342: $ ALPHA, BETA, VL, LDVL, VR, LDVR, WORK, -1,
343: $ RWORK, IERR )
344: LWKOPT = MAX( LWKOPT, N+INT( WORK( 1 ) ) )
345: ELSE
346: CALL ZGGHD3( JOBVL, JOBVR, N, 1, N, A, LDA, B, LDB, VL,
347: $ LDVL, VR, LDVR, WORK, -1, IERR )
348: LWKOPT = MAX( LWKOPT, N+INT( WORK( 1 ) ) )
349: CALL ZHGEQZ( 'E', JOBVL, JOBVR, N, 1, N, A, LDA, B, LDB,
350: $ ALPHA, BETA, VL, LDVL, VR, LDVR, WORK, -1,
351: $ RWORK, IERR )
352: LWKOPT = MAX( LWKOPT, N+INT( WORK( 1 ) ) )
353: END IF
354: WORK( 1 ) = DCMPLX( LWKOPT )
355: END IF
356: *
357: IF( INFO.NE.0 ) THEN
358: CALL XERBLA( 'ZGGEV3 ', -INFO )
359: RETURN
360: ELSE IF( LQUERY ) THEN
361: RETURN
362: END IF
363: *
364: * Quick return if possible
365: *
366: IF( N.EQ.0 )
367: $ RETURN
368: *
369: * Get machine constants
370: *
371: EPS = DLAMCH( 'E' )*DLAMCH( 'B' )
372: SMLNUM = DLAMCH( 'S' )
373: BIGNUM = ONE / SMLNUM
374: CALL DLABAD( SMLNUM, BIGNUM )
375: SMLNUM = SQRT( SMLNUM ) / EPS
376: BIGNUM = ONE / SMLNUM
377: *
378: * Scale A if max element outside range [SMLNUM,BIGNUM]
379: *
380: ANRM = ZLANGE( 'M', N, N, A, LDA, RWORK )
381: ILASCL = .FALSE.
382: IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
383: ANRMTO = SMLNUM
384: ILASCL = .TRUE.
385: ELSE IF( ANRM.GT.BIGNUM ) THEN
386: ANRMTO = BIGNUM
387: ILASCL = .TRUE.
388: END IF
389: IF( ILASCL )
390: $ CALL ZLASCL( 'G', 0, 0, ANRM, ANRMTO, N, N, A, LDA, IERR )
391: *
392: * Scale B if max element outside range [SMLNUM,BIGNUM]
393: *
394: BNRM = ZLANGE( 'M', N, N, B, LDB, RWORK )
395: ILBSCL = .FALSE.
396: IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
397: BNRMTO = SMLNUM
398: ILBSCL = .TRUE.
399: ELSE IF( BNRM.GT.BIGNUM ) THEN
400: BNRMTO = BIGNUM
401: ILBSCL = .TRUE.
402: END IF
403: IF( ILBSCL )
404: $ CALL ZLASCL( 'G', 0, 0, BNRM, BNRMTO, N, N, B, LDB, IERR )
405: *
406: * Permute the matrices A, B to isolate eigenvalues if possible
407: *
408: ILEFT = 1
409: IRIGHT = N + 1
410: IRWRK = IRIGHT + N
411: CALL ZGGBAL( 'P', N, A, LDA, B, LDB, ILO, IHI, RWORK( ILEFT ),
412: $ RWORK( IRIGHT ), RWORK( IRWRK ), IERR )
413: *
414: * Reduce B to triangular form (QR decomposition of B)
415: *
416: IROWS = IHI + 1 - ILO
417: IF( ILV ) THEN
418: ICOLS = N + 1 - ILO
419: ELSE
420: ICOLS = IROWS
421: END IF
422: ITAU = 1
423: IWRK = ITAU + IROWS
424: CALL ZGEQRF( IROWS, ICOLS, B( ILO, ILO ), LDB, WORK( ITAU ),
425: $ WORK( IWRK ), LWORK+1-IWRK, IERR )
426: *
427: * Apply the orthogonal transformation to matrix A
428: *
429: CALL ZUNMQR( 'L', 'C', IROWS, ICOLS, IROWS, B( ILO, ILO ), LDB,
430: $ WORK( ITAU ), A( ILO, ILO ), LDA, WORK( IWRK ),
431: $ LWORK+1-IWRK, IERR )
432: *
433: * Initialize VL
434: *
435: IF( ILVL ) THEN
436: CALL ZLASET( 'Full', N, N, CZERO, CONE, VL, LDVL )
437: IF( IROWS.GT.1 ) THEN
438: CALL ZLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB,
439: $ VL( ILO+1, ILO ), LDVL )
440: END IF
441: CALL ZUNGQR( IROWS, IROWS, IROWS, VL( ILO, ILO ), LDVL,
442: $ WORK( ITAU ), WORK( IWRK ), LWORK+1-IWRK, IERR )
443: END IF
444: *
445: * Initialize VR
446: *
447: IF( ILVR )
448: $ CALL ZLASET( 'Full', N, N, CZERO, CONE, VR, LDVR )
449: *
450: * Reduce to generalized Hessenberg form
451: *
452: IF( ILV ) THEN
453: *
454: * Eigenvectors requested -- work on whole matrix.
455: *
456: CALL ZGGHD3( JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB, VL,
457: $ LDVL, VR, LDVR, WORK( IWRK ), LWORK+1-IWRK, IERR )
458: ELSE
459: CALL ZGGHD3( 'N', 'N', IROWS, 1, IROWS, A( ILO, ILO ), LDA,
460: $ B( ILO, ILO ), LDB, VL, LDVL, VR, LDVR,
461: $ WORK( IWRK ), LWORK+1-IWRK, IERR )
462: END IF
463: *
464: * Perform QZ algorithm (Compute eigenvalues, and optionally, the
465: * Schur form and Schur vectors)
466: *
467: IWRK = ITAU
468: IF( ILV ) THEN
469: CHTEMP = 'S'
470: ELSE
471: CHTEMP = 'E'
472: END IF
473: CALL ZHGEQZ( CHTEMP, JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB,
474: $ ALPHA, BETA, VL, LDVL, VR, LDVR, WORK( IWRK ),
475: $ LWORK+1-IWRK, RWORK( IRWRK ), IERR )
476: IF( IERR.NE.0 ) THEN
477: IF( IERR.GT.0 .AND. IERR.LE.N ) THEN
478: INFO = IERR
479: ELSE IF( IERR.GT.N .AND. IERR.LE.2*N ) THEN
480: INFO = IERR - N
481: ELSE
482: INFO = N + 1
483: END IF
484: GO TO 70
485: END IF
486: *
487: * Compute Eigenvectors
488: *
489: IF( ILV ) THEN
490: IF( ILVL ) THEN
491: IF( ILVR ) THEN
492: CHTEMP = 'B'
493: ELSE
494: CHTEMP = 'L'
495: END IF
496: ELSE
497: CHTEMP = 'R'
498: END IF
499: *
500: CALL ZTGEVC( CHTEMP, 'B', LDUMMA, N, A, LDA, B, LDB, VL, LDVL,
501: $ VR, LDVR, N, IN, WORK( IWRK ), RWORK( IRWRK ),
502: $ IERR )
503: IF( IERR.NE.0 ) THEN
504: INFO = N + 2
505: GO TO 70
506: END IF
507: *
508: * Undo balancing on VL and VR and normalization
509: *
510: IF( ILVL ) THEN
511: CALL ZGGBAK( 'P', 'L', N, ILO, IHI, RWORK( ILEFT ),
512: $ RWORK( IRIGHT ), N, VL, LDVL, IERR )
513: DO 30 JC = 1, N
514: TEMP = ZERO
515: DO 10 JR = 1, N
516: TEMP = MAX( TEMP, ABS1( VL( JR, JC ) ) )
517: 10 CONTINUE
518: IF( TEMP.LT.SMLNUM )
519: $ GO TO 30
520: TEMP = ONE / TEMP
521: DO 20 JR = 1, N
522: VL( JR, JC ) = VL( JR, JC )*TEMP
523: 20 CONTINUE
524: 30 CONTINUE
525: END IF
526: IF( ILVR ) THEN
527: CALL ZGGBAK( 'P', 'R', N, ILO, IHI, RWORK( ILEFT ),
528: $ RWORK( IRIGHT ), N, VR, LDVR, IERR )
529: DO 60 JC = 1, N
530: TEMP = ZERO
531: DO 40 JR = 1, N
532: TEMP = MAX( TEMP, ABS1( VR( JR, JC ) ) )
533: 40 CONTINUE
534: IF( TEMP.LT.SMLNUM )
535: $ GO TO 60
536: TEMP = ONE / TEMP
537: DO 50 JR = 1, N
538: VR( JR, JC ) = VR( JR, JC )*TEMP
539: 50 CONTINUE
540: 60 CONTINUE
541: END IF
542: END IF
543: *
544: * Undo scaling if necessary
545: *
546: 70 CONTINUE
547: *
548: IF( ILASCL )
549: $ CALL ZLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHA, N, IERR )
550: *
551: IF( ILBSCL )
552: $ CALL ZLASCL( 'G', 0, 0, BNRMTO, BNRM, N, 1, BETA, N, IERR )
553: *
554: WORK( 1 ) = DCMPLX( LWKOPT )
555: RETURN
556: *
557: * End of ZGGEV3
558: *
559: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zggev3.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack