14 typedef long long BLASLONG;
15 typedef unsigned long long BLASULONG;
17 typedef long BLASLONG;
18 typedef unsigned long BLASULONG;
22 typedef BLASLONG blasint;
24 #define blasabs(x) llabs(x)
26 #define blasabs(x) labs(x)
30 #define blasabs(x) abs(x)
33 typedef blasint integer;
35 typedef unsigned int uinteger;
36 typedef char *address;
37 typedef short int shortint;
39 typedef double doublereal;
40 typedef struct { real r, i; } complex;
41 typedef struct { doublereal r, i; } doublecomplex;
43 static inline _Fcomplex Cf(complex *z) {_Fcomplex zz={z->r , z->i}; return zz;}
44 static inline _Dcomplex Cd(doublecomplex *z) {_Dcomplex zz={z->r , z->i};return zz;}
45 static inline _Fcomplex * _pCf(complex *z) {return (_Fcomplex*)z;}
46 static inline _Dcomplex * _pCd(doublecomplex *z) {return (_Dcomplex*)z;}
48 static inline _Complex float Cf(complex *z) {return z->r + z->i*_Complex_I;}
49 static inline _Complex double Cd(doublecomplex *z) {return z->r + z->i*_Complex_I;}
50 static inline _Complex float * _pCf(complex *z) {return (_Complex float*)z;}
51 static inline _Complex double * _pCd(doublecomplex *z) {return (_Complex double*)z;}
53 #define pCf(z) (*_pCf(z))
54 #define pCd(z) (*_pCd(z))
56 typedef short int shortlogical;
57 typedef char logical1;
58 typedef char integer1;
63 /* Extern is for use with -E */
74 /*external read, write*/
83 /*internal read, write*/
113 /*rewind, backspace, endfile*/
125 ftnint *inex; /*parameters in standard's order*/
151 union Multitype { /* for multiple entry points */
162 typedef union Multitype Multitype;
164 struct Vardesc { /* for Namelist */
170 typedef struct Vardesc Vardesc;
177 typedef struct Namelist Namelist;
179 #define abs(x) ((x) >= 0 ? (x) : -(x))
180 #define dabs(x) (fabs(x))
181 #define f2cmin(a,b) ((a) <= (b) ? (a) : (b))
182 #define f2cmax(a,b) ((a) >= (b) ? (a) : (b))
183 #define dmin(a,b) (f2cmin(a,b))
184 #define dmax(a,b) (f2cmax(a,b))
185 #define bit_test(a,b) ((a) >> (b) & 1)
186 #define bit_clear(a,b) ((a) & ~((uinteger)1 << (b)))
187 #define bit_set(a,b) ((a) | ((uinteger)1 << (b)))
189 #define abort_() { sig_die("Fortran abort routine called", 1); }
190 #define c_abs(z) (cabsf(Cf(z)))
191 #define c_cos(R,Z) { pCf(R)=ccos(Cf(Z)); }
193 #define c_div(c, a, b) {Cf(c)._Val[0] = (Cf(a)._Val[0]/Cf(b)._Val[0]); Cf(c)._Val[1]=(Cf(a)._Val[1]/Cf(b)._Val[1]);}
194 #define z_div(c, a, b) {Cd(c)._Val[0] = (Cd(a)._Val[0]/Cd(b)._Val[0]); Cd(c)._Val[1]=(Cd(a)._Val[1]/Cd(b)._Val[1]);}
196 #define c_div(c, a, b) {pCf(c) = Cf(a)/Cf(b);}
197 #define z_div(c, a, b) {pCd(c) = Cd(a)/Cd(b);}
199 #define c_exp(R, Z) {pCf(R) = cexpf(Cf(Z));}
200 #define c_log(R, Z) {pCf(R) = clogf(Cf(Z));}
201 #define c_sin(R, Z) {pCf(R) = csinf(Cf(Z));}
202 #define z_sin(R, Z) {pCd(R) = csin(Cd(Z));}
203 //#define c_sqrt(R, Z) {*(R) = csqrtf(Cf(Z));}
204 #define c_sqrt(R, Z) {pCf(R) = csqrtf(Cf(Z));}
205 #define d_abs(x) (fabs(*(x)))
206 #define d_acos(x) (acos(*(x)))
207 #define d_asin(x) (asin(*(x)))
208 #define d_atan(x) (atan(*(x)))
209 #define d_atn2(x, y) (atan2(*(x),*(y)))
210 #define d_cnjg(R, Z) { pCd(R) = conj(Cd(Z)); }
211 #define r_cnjg(R, Z) { pCf(R) = conjf(Cf(Z)); }
212 #define d_cos(x) (cos(*(x)))
213 #define d_cosh(x) (cosh(*(x)))
214 #define d_dim(__a, __b) ( *(__a) > *(__b) ? *(__a) - *(__b) : 0.0 )
215 #define d_exp(x) (exp(*(x)))
216 #define d_imag(z) (cimag(Cd(z)))
217 #define r_imag(z) (cimagf(Cf(z)))
218 #define d_int(__x) (*(__x)>0 ? floor(*(__x)) : -floor(- *(__x)))
219 #define r_int(__x) (*(__x)>0 ? floor(*(__x)) : -floor(- *(__x)))
220 #define d_lg10(x) ( 0.43429448190325182765 * log(*(x)) )
221 #define r_lg10(x) ( 0.43429448190325182765 * log(*(x)) )
222 #define d_log(x) (log(*(x)))
223 #define d_mod(x, y) (fmod(*(x), *(y)))
224 #define u_nint(__x) ((__x)>=0 ? floor((__x) + .5) : -floor(.5 - (__x)))
225 #define d_nint(x) u_nint(*(x))
226 #define u_sign(__a,__b) ((__b) >= 0 ? ((__a) >= 0 ? (__a) : -(__a)) : -((__a) >= 0 ? (__a) : -(__a)))
227 #define d_sign(a,b) u_sign(*(a),*(b))
228 #define r_sign(a,b) u_sign(*(a),*(b))
229 #define d_sin(x) (sin(*(x)))
230 #define d_sinh(x) (sinh(*(x)))
231 #define d_sqrt(x) (sqrt(*(x)))
232 #define d_tan(x) (tan(*(x)))
233 #define d_tanh(x) (tanh(*(x)))
234 #define i_abs(x) abs(*(x))
235 #define i_dnnt(x) ((integer)u_nint(*(x)))
236 #define i_len(s, n) (n)
237 #define i_nint(x) ((integer)u_nint(*(x)))
238 #define i_sign(a,b) ((integer)u_sign((integer)*(a),(integer)*(b)))
239 #define pow_dd(ap, bp) ( pow(*(ap), *(bp)))
240 #define pow_si(B,E) spow_ui(*(B),*(E))
241 #define pow_ri(B,E) spow_ui(*(B),*(E))
242 #define pow_di(B,E) dpow_ui(*(B),*(E))
243 #define pow_zi(p, a, b) {pCd(p) = zpow_ui(Cd(a), *(b));}
244 #define pow_ci(p, a, b) {pCf(p) = cpow_ui(Cf(a), *(b));}
245 #define pow_zz(R,A,B) {pCd(R) = cpow(Cd(A),*(B));}
246 #define s_cat(lpp, rpp, rnp, np, llp) { ftnlen i, nc, ll; char *f__rp, *lp; ll = (llp); lp = (lpp); for(i=0; i < (int)*(np); ++i) { nc = ll; if((rnp)[i] < nc) nc = (rnp)[i]; ll -= nc; f__rp = (rpp)[i]; while(--nc >= 0) *lp++ = *(f__rp)++; } while(--ll >= 0) *lp++ = ' '; }
247 #define s_cmp(a,b,c,d) ((integer)strncmp((a),(b),f2cmin((c),(d))))
248 #define s_copy(A,B,C,D) { int __i,__m; for (__i=0, __m=f2cmin((C),(D)); __i<__m && (B)[__i] != 0; ++__i) (A)[__i] = (B)[__i]; }
249 #define sig_die(s, kill) { exit(1); }
250 #define s_stop(s, n) {exit(0);}
251 static char junk[] = "\n@(#)LIBF77 VERSION 19990503\n";
252 #define z_abs(z) (cabs(Cd(z)))
253 #define z_exp(R, Z) {pCd(R) = cexp(Cd(Z));}
254 #define z_sqrt(R, Z) {pCd(R) = csqrt(Cd(Z));}
255 #define myexit_() break;
256 #define mycycle_() continue;
257 #define myceiling_(w) {ceil(w)}
258 #define myhuge_(w) {HUGE_VAL}
259 //#define mymaxloc_(w,s,e,n) {if (sizeof(*(w)) == sizeof(double)) dmaxloc_((w),*(s),*(e),n); else dmaxloc_((w),*(s),*(e),n);}
260 #define mymaxloc_(w,s,e,n) {dmaxloc_(w,*(s),*(e),n)}
262 /* procedure parameter types for -A and -C++ */
264 #define F2C_proc_par_types 1
266 typedef logical (*L_fp)(...);
268 typedef logical (*L_fp)();
271 static float spow_ui(float x, integer n) {
272 float pow=1.0; unsigned long int u;
274 if(n < 0) n = -n, x = 1/x;
283 static double dpow_ui(double x, integer n) {
284 double pow=1.0; unsigned long int u;
286 if(n < 0) n = -n, x = 1/x;
296 static _Fcomplex cpow_ui(complex x, integer n) {
297 complex pow={1.0,0.0}; unsigned long int u;
299 if(n < 0) n = -n, x.r = 1/x.r, x.i=1/x.i;
301 if(u & 01) pow.r *= x.r, pow.i *= x.i;
302 if(u >>= 1) x.r *= x.r, x.i *= x.i;
306 _Fcomplex p={pow.r, pow.i};
310 static _Complex float cpow_ui(_Complex float x, integer n) {
311 _Complex float pow=1.0; unsigned long int u;
313 if(n < 0) n = -n, x = 1/x;
324 static _Dcomplex zpow_ui(_Dcomplex x, integer n) {
325 _Dcomplex pow={1.0,0.0}; unsigned long int u;
327 if(n < 0) n = -n, x._Val[0] = 1/x._Val[0], x._Val[1] =1/x._Val[1];
329 if(u & 01) pow._Val[0] *= x._Val[0], pow._Val[1] *= x._Val[1];
330 if(u >>= 1) x._Val[0] *= x._Val[0], x._Val[1] *= x._Val[1];
334 _Dcomplex p = {pow._Val[0], pow._Val[1]};
338 static _Complex double zpow_ui(_Complex double x, integer n) {
339 _Complex double pow=1.0; unsigned long int u;
341 if(n < 0) n = -n, x = 1/x;
351 static integer pow_ii(integer x, integer n) {
352 integer pow; unsigned long int u;
354 if (n == 0 || x == 1) pow = 1;
355 else if (x != -1) pow = x == 0 ? 1/x : 0;
358 if ((n > 0) || !(n == 0 || x == 1 || x != -1)) {
368 static integer dmaxloc_(double *w, integer s, integer e, integer *n)
370 double m; integer i, mi;
371 for(m=w[s-1], mi=s, i=s+1; i<=e; i++)
372 if (w[i-1]>m) mi=i ,m=w[i-1];
375 static integer smaxloc_(float *w, integer s, integer e, integer *n)
377 float m; integer i, mi;
378 for(m=w[s-1], mi=s, i=s+1; i<=e; i++)
379 if (w[i-1]>m) mi=i ,m=w[i-1];
382 static inline void cdotc_(complex *z, integer *n_, complex *x, integer *incx_, complex *y, integer *incy_) {
383 integer n = *n_, incx = *incx_, incy = *incy_, i;
385 _Fcomplex zdotc = {0.0, 0.0};
386 if (incx == 1 && incy == 1) {
387 for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
388 zdotc._Val[0] += conjf(Cf(&x[i]))._Val[0] * Cf(&y[i])._Val[0];
389 zdotc._Val[1] += conjf(Cf(&x[i]))._Val[1] * Cf(&y[i])._Val[1];
392 for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
393 zdotc._Val[0] += conjf(Cf(&x[i*incx]))._Val[0] * Cf(&y[i*incy])._Val[0];
394 zdotc._Val[1] += conjf(Cf(&x[i*incx]))._Val[1] * Cf(&y[i*incy])._Val[1];
400 _Complex float zdotc = 0.0;
401 if (incx == 1 && incy == 1) {
402 for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
403 zdotc += conjf(Cf(&x[i])) * Cf(&y[i]);
406 for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
407 zdotc += conjf(Cf(&x[i*incx])) * Cf(&y[i*incy]);
413 static inline void zdotc_(doublecomplex *z, integer *n_, doublecomplex *x, integer *incx_, doublecomplex *y, integer *incy_) {
414 integer n = *n_, incx = *incx_, incy = *incy_, i;
416 _Dcomplex zdotc = {0.0, 0.0};
417 if (incx == 1 && incy == 1) {
418 for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
419 zdotc._Val[0] += conj(Cd(&x[i]))._Val[0] * Cd(&y[i])._Val[0];
420 zdotc._Val[1] += conj(Cd(&x[i]))._Val[1] * Cd(&y[i])._Val[1];
423 for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
424 zdotc._Val[0] += conj(Cd(&x[i*incx]))._Val[0] * Cd(&y[i*incy])._Val[0];
425 zdotc._Val[1] += conj(Cd(&x[i*incx]))._Val[1] * Cd(&y[i*incy])._Val[1];
431 _Complex double zdotc = 0.0;
432 if (incx == 1 && incy == 1) {
433 for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
434 zdotc += conj(Cd(&x[i])) * Cd(&y[i]);
437 for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
438 zdotc += conj(Cd(&x[i*incx])) * Cd(&y[i*incy]);
444 static inline void cdotu_(complex *z, integer *n_, complex *x, integer *incx_, complex *y, integer *incy_) {
445 integer n = *n_, incx = *incx_, incy = *incy_, i;
447 _Fcomplex zdotc = {0.0, 0.0};
448 if (incx == 1 && incy == 1) {
449 for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
450 zdotc._Val[0] += Cf(&x[i])._Val[0] * Cf(&y[i])._Val[0];
451 zdotc._Val[1] += Cf(&x[i])._Val[1] * Cf(&y[i])._Val[1];
454 for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
455 zdotc._Val[0] += Cf(&x[i*incx])._Val[0] * Cf(&y[i*incy])._Val[0];
456 zdotc._Val[1] += Cf(&x[i*incx])._Val[1] * Cf(&y[i*incy])._Val[1];
462 _Complex float zdotc = 0.0;
463 if (incx == 1 && incy == 1) {
464 for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
465 zdotc += Cf(&x[i]) * Cf(&y[i]);
468 for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
469 zdotc += Cf(&x[i*incx]) * Cf(&y[i*incy]);
475 static inline void zdotu_(doublecomplex *z, integer *n_, doublecomplex *x, integer *incx_, doublecomplex *y, integer *incy_) {
476 integer n = *n_, incx = *incx_, incy = *incy_, i;
478 _Dcomplex zdotc = {0.0, 0.0};
479 if (incx == 1 && incy == 1) {
480 for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
481 zdotc._Val[0] += Cd(&x[i])._Val[0] * Cd(&y[i])._Val[0];
482 zdotc._Val[1] += Cd(&x[i])._Val[1] * Cd(&y[i])._Val[1];
485 for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
486 zdotc._Val[0] += Cd(&x[i*incx])._Val[0] * Cd(&y[i*incy])._Val[0];
487 zdotc._Val[1] += Cd(&x[i*incx])._Val[1] * Cd(&y[i*incy])._Val[1];
493 _Complex double zdotc = 0.0;
494 if (incx == 1 && incy == 1) {
495 for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
496 zdotc += Cd(&x[i]) * Cd(&y[i]);
499 for (i=0;i<n;i++) { /* zdotc = zdotc + dconjg(x(i))* y(i) */
500 zdotc += Cd(&x[i*incx]) * Cd(&y[i*incy]);
506 /* -- translated by f2c (version 20000121).
507 You must link the resulting object file with the libraries:
508 -lf2c -lm (in that order)
515 /* Table of constant values */
517 static doublecomplex c_b1 = {1.,0.};
518 static doublecomplex c_b3 = {0.,0.};
519 static doublecomplex c_b5 = {20.,0.};
521 /* > \brief \b ZLATM5 */
523 /* =========== DOCUMENTATION =========== */
525 /* Online html documentation available at */
526 /* http://www.netlib.org/lapack/explore-html/ */
531 /* SUBROUTINE ZLATM5( PRTYPE, M, N, A, LDA, B, LDB, C, LDC, D, LDD, */
532 /* E, LDE, F, LDF, R, LDR, L, LDL, ALPHA, QBLCKA, */
535 /* INTEGER LDA, LDB, LDC, LDD, LDE, LDF, LDL, LDR, M, N, */
536 /* $ PRTYPE, QBLCKA, QBLCKB */
537 /* DOUBLE PRECISION ALPHA */
538 /* COMPLEX*16 A( LDA, * ), B( LDB, * ), C( LDC, * ), */
539 /* $ D( LDD, * ), E( LDE, * ), F( LDF, * ), */
540 /* $ L( LDL, * ), R( LDR, * ) */
543 /* > \par Purpose: */
548 /* > ZLATM5 generates matrices involved in the Generalized Sylvester */
551 /* > A * R - L * B = C */
552 /* > D * R - L * E = F */
554 /* > They also satisfy (the diagonalization condition) */
556 /* > [ I -L ] ( [ A -C ], [ D -F ] ) [ I R ] = ( [ A ], [ D ] ) */
557 /* > [ I ] ( [ B ] [ E ] ) [ I ] ( [ B ] [ E ] ) */
564 /* > \param[in] PRTYPE */
566 /* > PRTYPE is INTEGER */
567 /* > "Points" to a certain type of the matrices to generate */
568 /* > (see further details). */
574 /* > Specifies the order of A and D and the number of rows in */
575 /* > C, F, R and L. */
581 /* > Specifies the order of B and E and the number of columns in */
582 /* > C, F, R and L. */
585 /* > \param[out] A */
587 /* > A is COMPLEX*16 array, dimension (LDA, M). */
588 /* > On exit A M-by-M is initialized according to PRTYPE. */
591 /* > \param[in] LDA */
593 /* > LDA is INTEGER */
594 /* > The leading dimension of A. */
597 /* > \param[out] B */
599 /* > B is COMPLEX*16 array, dimension (LDB, N). */
600 /* > On exit B N-by-N is initialized according to PRTYPE. */
603 /* > \param[in] LDB */
605 /* > LDB is INTEGER */
606 /* > The leading dimension of B. */
609 /* > \param[out] C */
611 /* > C is COMPLEX*16 array, dimension (LDC, N). */
612 /* > On exit C M-by-N is initialized according to PRTYPE. */
615 /* > \param[in] LDC */
617 /* > LDC is INTEGER */
618 /* > The leading dimension of C. */
621 /* > \param[out] D */
623 /* > D is COMPLEX*16 array, dimension (LDD, M). */
624 /* > On exit D M-by-M is initialized according to PRTYPE. */
627 /* > \param[in] LDD */
629 /* > LDD is INTEGER */
630 /* > The leading dimension of D. */
633 /* > \param[out] E */
635 /* > E is COMPLEX*16 array, dimension (LDE, N). */
636 /* > On exit E N-by-N is initialized according to PRTYPE. */
639 /* > \param[in] LDE */
641 /* > LDE is INTEGER */
642 /* > The leading dimension of E. */
645 /* > \param[out] F */
647 /* > F is COMPLEX*16 array, dimension (LDF, N). */
648 /* > On exit F M-by-N is initialized according to PRTYPE. */
651 /* > \param[in] LDF */
653 /* > LDF is INTEGER */
654 /* > The leading dimension of F. */
657 /* > \param[out] R */
659 /* > R is COMPLEX*16 array, dimension (LDR, N). */
660 /* > On exit R M-by-N is initialized according to PRTYPE. */
663 /* > \param[in] LDR */
665 /* > LDR is INTEGER */
666 /* > The leading dimension of R. */
669 /* > \param[out] L */
671 /* > L is COMPLEX*16 array, dimension (LDL, N). */
672 /* > On exit L M-by-N is initialized according to PRTYPE. */
675 /* > \param[in] LDL */
677 /* > LDL is INTEGER */
678 /* > The leading dimension of L. */
681 /* > \param[in] ALPHA */
683 /* > ALPHA is DOUBLE PRECISION */
684 /* > Parameter used in generating PRTYPE = 1 and 5 matrices. */
687 /* > \param[in] QBLCKA */
689 /* > QBLCKA is INTEGER */
690 /* > When PRTYPE = 3, specifies the distance between 2-by-2 */
691 /* > blocks on the diagonal in A. Otherwise, QBLCKA is not */
692 /* > referenced. QBLCKA > 1. */
695 /* > \param[in] QBLCKB */
697 /* > QBLCKB is INTEGER */
698 /* > When PRTYPE = 3, specifies the distance between 2-by-2 */
699 /* > blocks on the diagonal in B. Otherwise, QBLCKB is not */
700 /* > referenced. QBLCKB > 1. */
706 /* > \author Univ. of Tennessee */
707 /* > \author Univ. of California Berkeley */
708 /* > \author Univ. of Colorado Denver */
709 /* > \author NAG Ltd. */
711 /* > \date June 2016 */
713 /* > \ingroup complex16_matgen */
715 /* > \par Further Details: */
716 /* ===================== */
720 /* > PRTYPE = 1: A and B are Jordan blocks, D and E are identity matrices */
722 /* > A : if (i == j) then A(i, j) = 1.0 */
723 /* > if (j == i + 1) then A(i, j) = -1.0 */
724 /* > else A(i, j) = 0.0, i, j = 1...M */
726 /* > B : if (i == j) then B(i, j) = 1.0 - ALPHA */
727 /* > if (j == i + 1) then B(i, j) = 1.0 */
728 /* > else B(i, j) = 0.0, i, j = 1...N */
730 /* > D : if (i == j) then D(i, j) = 1.0 */
731 /* > else D(i, j) = 0.0, i, j = 1...M */
733 /* > E : if (i == j) then E(i, j) = 1.0 */
734 /* > else E(i, j) = 0.0, i, j = 1...N */
736 /* > L = R are chosen from [-10...10], */
737 /* > which specifies the right hand sides (C, F). */
739 /* > PRTYPE = 2 or 3: Triangular and/or quasi- triangular. */
741 /* > A : if (i <= j) then A(i, j) = [-1...1] */
742 /* > else A(i, j) = 0.0, i, j = 1...M */
744 /* > if (PRTYPE = 3) then */
745 /* > A(k + 1, k + 1) = A(k, k) */
746 /* > A(k + 1, k) = [-1...1] */
747 /* > sign(A(k, k + 1) = -(sin(A(k + 1, k)) */
748 /* > k = 1, M - 1, QBLCKA */
750 /* > B : if (i <= j) then B(i, j) = [-1...1] */
751 /* > else B(i, j) = 0.0, i, j = 1...N */
753 /* > if (PRTYPE = 3) then */
754 /* > B(k + 1, k + 1) = B(k, k) */
755 /* > B(k + 1, k) = [-1...1] */
756 /* > sign(B(k, k + 1) = -(sign(B(k + 1, k)) */
757 /* > k = 1, N - 1, QBLCKB */
759 /* > D : if (i <= j) then D(i, j) = [-1...1]. */
760 /* > else D(i, j) = 0.0, i, j = 1...M */
763 /* > E : if (i <= j) then D(i, j) = [-1...1] */
764 /* > else E(i, j) = 0.0, i, j = 1...N */
766 /* > L, R are chosen from [-10...10], */
767 /* > which specifies the right hand sides (C, F). */
769 /* > PRTYPE = 4 Full */
770 /* > A(i, j) = [-10...10] */
771 /* > D(i, j) = [-1...1] i,j = 1...M */
772 /* > B(i, j) = [-10...10] */
773 /* > E(i, j) = [-1...1] i,j = 1...N */
774 /* > R(i, j) = [-10...10] */
775 /* > L(i, j) = [-1...1] i = 1..M ,j = 1...N */
777 /* > L, R specifies the right hand sides (C, F). */
779 /* > PRTYPE = 5 special case common and/or close eigs. */
782 /* ===================================================================== */
783 /* Subroutine */ int zlatm5_(integer *prtype, integer *m, integer *n,
784 doublecomplex *a, integer *lda, doublecomplex *b, integer *ldb,
785 doublecomplex *c__, integer *ldc, doublecomplex *d__, integer *ldd,
786 doublecomplex *e, integer *lde, doublecomplex *f, integer *ldf,
787 doublecomplex *r__, integer *ldr, doublecomplex *l, integer *ldl,
788 doublereal *alpha, integer *qblcka, integer *qblckb)
790 /* System generated locals */
791 integer a_dim1, a_offset, b_dim1, b_offset, c_dim1, c_offset, d_dim1,
792 d_offset, e_dim1, e_offset, f_dim1, f_offset, l_dim1, l_offset,
793 r_dim1, r_offset, i__1, i__2, i__3, i__4;
795 doublecomplex z__1, z__2, z__3, z__4, z__5;
797 /* Local variables */
799 doublecomplex imeps, reeps;
800 extern /* Subroutine */ int zgemm_(char *, char *, integer *, integer *,
801 integer *, doublecomplex *, doublecomplex *, integer *,
802 doublecomplex *, integer *, doublecomplex *, doublecomplex *,
806 /* -- LAPACK computational routine (version 3.7.0) -- */
807 /* -- LAPACK is a software package provided by Univ. of Tennessee, -- */
808 /* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */
812 /* ===================================================================== */
815 /* Parameter adjustments */
817 a_offset = 1 + a_dim1 * 1;
820 b_offset = 1 + b_dim1 * 1;
823 c_offset = 1 + c_dim1 * 1;
826 d_offset = 1 + d_dim1 * 1;
829 e_offset = 1 + e_dim1 * 1;
832 f_offset = 1 + f_dim1 * 1;
835 r_offset = 1 + r_dim1 * 1;
838 l_offset = 1 + l_dim1 * 1;
844 for (i__ = 1; i__ <= i__1; ++i__) {
846 for (j = 1; j <= i__2; ++j) {
848 i__3 = i__ + j * a_dim1;
849 a[i__3].r = 1., a[i__3].i = 0.;
850 i__3 = i__ + j * d_dim1;
851 d__[i__3].r = 1., d__[i__3].i = 0.;
852 } else if (i__ == j - 1) {
853 i__3 = i__ + j * a_dim1;
854 z__1.r = -1., z__1.i = 0.;
855 a[i__3].r = z__1.r, a[i__3].i = z__1.i;
856 i__3 = i__ + j * d_dim1;
857 d__[i__3].r = 0., d__[i__3].i = 0.;
859 i__3 = i__ + j * a_dim1;
860 a[i__3].r = 0., a[i__3].i = 0.;
861 i__3 = i__ + j * d_dim1;
862 d__[i__3].r = 0., d__[i__3].i = 0.;
870 for (i__ = 1; i__ <= i__1; ++i__) {
872 for (j = 1; j <= i__2; ++j) {
874 i__3 = i__ + j * b_dim1;
875 z__1.r = 1. - *alpha, z__1.i = 0.;
876 b[i__3].r = z__1.r, b[i__3].i = z__1.i;
877 i__3 = i__ + j * e_dim1;
878 e[i__3].r = 1., e[i__3].i = 0.;
879 } else if (i__ == j - 1) {
880 i__3 = i__ + j * b_dim1;
881 b[i__3].r = 1., b[i__3].i = 0.;
882 i__3 = i__ + j * e_dim1;
883 e[i__3].r = 0., e[i__3].i = 0.;
885 i__3 = i__ + j * b_dim1;
886 b[i__3].r = 0., b[i__3].i = 0.;
887 i__3 = i__ + j * e_dim1;
888 e[i__3].r = 0., e[i__3].i = 0.;
896 for (i__ = 1; i__ <= i__1; ++i__) {
898 for (j = 1; j <= i__2; ++j) {
899 i__3 = i__ + j * r_dim1;
901 z__4.r = (doublereal) i__4, z__4.i = 0.;
903 z__2.r = .5 - z__3.r, z__2.i = 0. - z__3.i;
904 z__1.r = z__2.r * 20. - z__2.i * 0., z__1.i = z__2.r * 0. +
906 r__[i__3].r = z__1.r, r__[i__3].i = z__1.i;
907 i__3 = i__ + j * l_dim1;
908 i__4 = i__ + j * r_dim1;
909 l[i__3].r = r__[i__4].r, l[i__3].i = r__[i__4].i;
915 } else if (*prtype == 2 || *prtype == 3) {
917 for (i__ = 1; i__ <= i__1; ++i__) {
919 for (j = 1; j <= i__2; ++j) {
921 i__3 = i__ + j * a_dim1;
922 z__4.r = (doublereal) i__, z__4.i = 0.;
924 z__2.r = .5 - z__3.r, z__2.i = 0. - z__3.i;
925 z__1.r = z__2.r * 2. - z__2.i * 0., z__1.i = z__2.r * 0.
927 a[i__3].r = z__1.r, a[i__3].i = z__1.i;
928 i__3 = i__ + j * d_dim1;
930 z__4.r = (doublereal) i__4, z__4.i = 0.;
932 z__2.r = .5 - z__3.r, z__2.i = 0. - z__3.i;
933 z__1.r = z__2.r * 2. - z__2.i * 0., z__1.i = z__2.r * 0.
935 d__[i__3].r = z__1.r, d__[i__3].i = z__1.i;
937 i__3 = i__ + j * a_dim1;
938 a[i__3].r = 0., a[i__3].i = 0.;
939 i__3 = i__ + j * d_dim1;
940 d__[i__3].r = 0., d__[i__3].i = 0.;
948 for (i__ = 1; i__ <= i__1; ++i__) {
950 for (j = 1; j <= i__2; ++j) {
952 i__3 = i__ + j * b_dim1;
954 z__4.r = (doublereal) i__4, z__4.i = 0.;
956 z__2.r = .5 - z__3.r, z__2.i = 0. - z__3.i;
957 z__1.r = z__2.r * 2. - z__2.i * 0., z__1.i = z__2.r * 0.
959 b[i__3].r = z__1.r, b[i__3].i = z__1.i;
960 i__3 = i__ + j * e_dim1;
961 z__4.r = (doublereal) j, z__4.i = 0.;
963 z__2.r = .5 - z__3.r, z__2.i = 0. - z__3.i;
964 z__1.r = z__2.r * 2. - z__2.i * 0., z__1.i = z__2.r * 0.
966 e[i__3].r = z__1.r, e[i__3].i = z__1.i;
968 i__3 = i__ + j * b_dim1;
969 b[i__3].r = 0., b[i__3].i = 0.;
970 i__3 = i__ + j * e_dim1;
971 e[i__3].r = 0., e[i__3].i = 0.;
979 for (i__ = 1; i__ <= i__1; ++i__) {
981 for (j = 1; j <= i__2; ++j) {
982 i__3 = i__ + j * r_dim1;
984 z__4.r = (doublereal) i__4, z__4.i = 0.;
986 z__2.r = .5 - z__3.r, z__2.i = 0. - z__3.i;
987 z__1.r = z__2.r * 20. - z__2.i * 0., z__1.i = z__2.r * 0. +
989 r__[i__3].r = z__1.r, r__[i__3].i = z__1.i;
990 i__3 = i__ + j * l_dim1;
992 z__4.r = (doublereal) i__4, z__4.i = 0.;
994 z__2.r = .5 - z__3.r, z__2.i = 0. - z__3.i;
995 z__1.r = z__2.r * 20. - z__2.i * 0., z__1.i = z__2.r * 0. +
997 l[i__3].r = z__1.r, l[i__3].i = z__1.i;
1009 for (k = 1; i__2 < 0 ? k >= i__1 : k <= i__1; k += i__2) {
1010 i__3 = k + 1 + (k + 1) * a_dim1;
1011 i__4 = k + k * a_dim1;
1012 a[i__3].r = a[i__4].r, a[i__3].i = a[i__4].i;
1013 i__3 = k + 1 + k * a_dim1;
1014 z_sin(&z__2, &a[k + (k + 1) * a_dim1]);
1015 z__1.r = -z__2.r, z__1.i = -z__2.i;
1016 a[i__3].r = z__1.r, a[i__3].i = z__1.i;
1025 for (k = 1; i__1 < 0 ? k >= i__2 : k <= i__2; k += i__1) {
1026 i__3 = k + 1 + (k + 1) * b_dim1;
1027 i__4 = k + k * b_dim1;
1028 b[i__3].r = b[i__4].r, b[i__3].i = b[i__4].i;
1029 i__3 = k + 1 + k * b_dim1;
1030 z_sin(&z__2, &b[k + (k + 1) * b_dim1]);
1031 z__1.r = -z__2.r, z__1.i = -z__2.i;
1032 b[i__3].r = z__1.r, b[i__3].i = z__1.i;
1037 } else if (*prtype == 4) {
1039 for (i__ = 1; i__ <= i__1; ++i__) {
1041 for (j = 1; j <= i__2; ++j) {
1042 i__3 = i__ + j * a_dim1;
1044 z__4.r = (doublereal) i__4, z__4.i = 0.;
1045 z_sin(&z__3, &z__4);
1046 z__2.r = .5 - z__3.r, z__2.i = 0. - z__3.i;
1047 z__1.r = z__2.r * 20. - z__2.i * 0., z__1.i = z__2.r * 0. +
1049 a[i__3].r = z__1.r, a[i__3].i = z__1.i;
1050 i__3 = i__ + j * d_dim1;
1052 z__4.r = (doublereal) i__4, z__4.i = 0.;
1053 z_sin(&z__3, &z__4);
1054 z__2.r = .5 - z__3.r, z__2.i = 0. - z__3.i;
1055 z__1.r = z__2.r * 2. - z__2.i * 0., z__1.i = z__2.r * 0. +
1057 d__[i__3].r = z__1.r, d__[i__3].i = z__1.i;
1064 for (i__ = 1; i__ <= i__1; ++i__) {
1066 for (j = 1; j <= i__2; ++j) {
1067 i__3 = i__ + j * b_dim1;
1069 z__4.r = (doublereal) i__4, z__4.i = 0.;
1070 z_sin(&z__3, &z__4);
1071 z__2.r = .5 - z__3.r, z__2.i = 0. - z__3.i;
1072 z__1.r = z__2.r * 20. - z__2.i * 0., z__1.i = z__2.r * 0. +
1074 b[i__3].r = z__1.r, b[i__3].i = z__1.i;
1075 i__3 = i__ + j * e_dim1;
1077 z__4.r = (doublereal) i__4, z__4.i = 0.;
1078 z_sin(&z__3, &z__4);
1079 z__2.r = .5 - z__3.r, z__2.i = 0. - z__3.i;
1080 z__1.r = z__2.r * 2. - z__2.i * 0., z__1.i = z__2.r * 0. +
1082 e[i__3].r = z__1.r, e[i__3].i = z__1.i;
1089 for (i__ = 1; i__ <= i__1; ++i__) {
1091 for (j = 1; j <= i__2; ++j) {
1092 i__3 = i__ + j * r_dim1;
1094 z__4.r = (doublereal) i__4, z__4.i = 0.;
1095 z_sin(&z__3, &z__4);
1096 z__2.r = .5 - z__3.r, z__2.i = 0. - z__3.i;
1097 z__1.r = z__2.r * 20. - z__2.i * 0., z__1.i = z__2.r * 0. +
1099 r__[i__3].r = z__1.r, r__[i__3].i = z__1.i;
1100 i__3 = i__ + j * l_dim1;
1102 z__4.r = (doublereal) i__4, z__4.i = 0.;
1103 z_sin(&z__3, &z__4);
1104 z__2.r = .5 - z__3.r, z__2.i = 0. - z__3.i;
1105 z__1.r = z__2.r * 2. - z__2.i * 0., z__1.i = z__2.r * 0. +
1107 l[i__3].r = z__1.r, l[i__3].i = z__1.i;
1113 } else if (*prtype >= 5) {
1114 z__3.r = 1., z__3.i = 0.;
1115 z__2.r = z__3.r * 20. - z__3.i * 0., z__2.i = z__3.r * 0. + z__3.i *
1117 z__1.r = z__2.r / *alpha, z__1.i = z__2.i / *alpha;
1118 reeps.r = z__1.r, reeps.i = z__1.i;
1119 z__2.r = -1.5, z__2.i = 0.;
1120 z__1.r = z__2.r / *alpha, z__1.i = z__2.i / *alpha;
1121 imeps.r = z__1.r, imeps.i = z__1.i;
1123 for (i__ = 1; i__ <= i__1; ++i__) {
1125 for (j = 1; j <= i__2; ++j) {
1126 i__3 = i__ + j * r_dim1;
1128 z__5.r = (doublereal) i__4, z__5.i = 0.;
1129 z_sin(&z__4, &z__5);
1130 z__3.r = .5 - z__4.r, z__3.i = 0. - z__4.i;
1131 z__2.r = *alpha * z__3.r, z__2.i = *alpha * z__3.i;
1132 z_div(&z__1, &z__2, &c_b5);
1133 r__[i__3].r = z__1.r, r__[i__3].i = z__1.i;
1134 i__3 = i__ + j * l_dim1;
1136 z__5.r = (doublereal) i__4, z__5.i = 0.;
1137 z_sin(&z__4, &z__5);
1138 z__3.r = .5 - z__4.r, z__3.i = 0. - z__4.i;
1139 z__2.r = *alpha * z__3.r, z__2.i = *alpha * z__3.i;
1140 z_div(&z__1, &z__2, &c_b5);
1141 l[i__3].r = z__1.r, l[i__3].i = z__1.i;
1148 for (i__ = 1; i__ <= i__1; ++i__) {
1149 i__2 = i__ + i__ * d_dim1;
1150 d__[i__2].r = 1., d__[i__2].i = 0.;
1155 for (i__ = 1; i__ <= i__1; ++i__) {
1157 i__2 = i__ + i__ * a_dim1;
1158 a[i__2].r = 1., a[i__2].i = 0.;
1160 i__2 = i__ + i__ * a_dim1;
1161 z__1.r = reeps.r + 1., z__1.i = reeps.i + 0.;
1162 a[i__2].r = z__1.r, a[i__2].i = z__1.i;
1164 if (i__ % 2 != 0 && i__ < *m) {
1165 i__2 = i__ + (i__ + 1) * a_dim1;
1166 a[i__2].r = imeps.r, a[i__2].i = imeps.i;
1167 } else if (i__ > 1) {
1168 i__2 = i__ + (i__ - 1) * a_dim1;
1169 z__1.r = -imeps.r, z__1.i = -imeps.i;
1170 a[i__2].r = z__1.r, a[i__2].i = z__1.i;
1172 } else if (i__ <= 8) {
1174 i__2 = i__ + i__ * a_dim1;
1175 a[i__2].r = reeps.r, a[i__2].i = reeps.i;
1177 i__2 = i__ + i__ * a_dim1;
1178 z__1.r = -reeps.r, z__1.i = -reeps.i;
1179 a[i__2].r = z__1.r, a[i__2].i = z__1.i;
1181 if (i__ % 2 != 0 && i__ < *m) {
1182 i__2 = i__ + (i__ + 1) * a_dim1;
1183 a[i__2].r = 1., a[i__2].i = 0.;
1184 } else if (i__ > 1) {
1185 i__2 = i__ + (i__ - 1) * a_dim1;
1186 z__1.r = -1., z__1.i = 0.;
1187 a[i__2].r = z__1.r, a[i__2].i = z__1.i;
1190 i__2 = i__ + i__ * a_dim1;
1191 a[i__2].r = 1., a[i__2].i = 0.;
1192 if (i__ % 2 != 0 && i__ < *m) {
1193 i__2 = i__ + (i__ + 1) * a_dim1;
1195 z__1.r = d__1 * imeps.r, z__1.i = d__1 * imeps.i;
1196 a[i__2].r = z__1.r, a[i__2].i = z__1.i;
1197 } else if (i__ > 1) {
1198 i__2 = i__ + (i__ - 1) * a_dim1;
1199 z__2.r = -imeps.r, z__2.i = -imeps.i;
1201 z__1.r = d__1 * z__2.r, z__1.i = d__1 * z__2.i;
1202 a[i__2].r = z__1.r, a[i__2].i = z__1.i;
1209 for (i__ = 1; i__ <= i__1; ++i__) {
1210 i__2 = i__ + i__ * e_dim1;
1211 e[i__2].r = 1., e[i__2].i = 0.;
1213 i__2 = i__ + i__ * b_dim1;
1214 z__1.r = -1., z__1.i = 0.;
1215 b[i__2].r = z__1.r, b[i__2].i = z__1.i;
1217 i__2 = i__ + i__ * b_dim1;
1218 z__1.r = 1. - reeps.r, z__1.i = 0. - reeps.i;
1219 b[i__2].r = z__1.r, b[i__2].i = z__1.i;
1221 if (i__ % 2 != 0 && i__ < *n) {
1222 i__2 = i__ + (i__ + 1) * b_dim1;
1223 b[i__2].r = imeps.r, b[i__2].i = imeps.i;
1224 } else if (i__ > 1) {
1225 i__2 = i__ + (i__ - 1) * b_dim1;
1226 z__1.r = -imeps.r, z__1.i = -imeps.i;
1227 b[i__2].r = z__1.r, b[i__2].i = z__1.i;
1229 } else if (i__ <= 8) {
1231 i__2 = i__ + i__ * b_dim1;
1232 b[i__2].r = reeps.r, b[i__2].i = reeps.i;
1234 i__2 = i__ + i__ * b_dim1;
1235 z__1.r = -reeps.r, z__1.i = -reeps.i;
1236 b[i__2].r = z__1.r, b[i__2].i = z__1.i;
1238 if (i__ % 2 != 0 && i__ < *n) {
1239 i__2 = i__ + (i__ + 1) * b_dim1;
1240 z__1.r = imeps.r + 1., z__1.i = imeps.i + 0.;
1241 b[i__2].r = z__1.r, b[i__2].i = z__1.i;
1242 } else if (i__ > 1) {
1243 i__2 = i__ + (i__ - 1) * b_dim1;
1244 z__2.r = -1., z__2.i = 0.;
1245 z__1.r = z__2.r - imeps.r, z__1.i = z__2.i - imeps.i;
1246 b[i__2].r = z__1.r, b[i__2].i = z__1.i;
1249 i__2 = i__ + i__ * b_dim1;
1250 z__1.r = 1. - reeps.r, z__1.i = 0. - reeps.i;
1251 b[i__2].r = z__1.r, b[i__2].i = z__1.i;
1252 if (i__ % 2 != 0 && i__ < *n) {
1253 i__2 = i__ + (i__ + 1) * b_dim1;
1255 z__1.r = d__1 * imeps.r, z__1.i = d__1 * imeps.i;
1256 b[i__2].r = z__1.r, b[i__2].i = z__1.i;
1257 } else if (i__ > 1) {
1258 i__2 = i__ + (i__ - 1) * b_dim1;
1259 z__2.r = -imeps.r, z__2.i = -imeps.i;
1261 z__1.r = d__1 * z__2.r, z__1.i = d__1 * z__2.i;
1262 b[i__2].r = z__1.r, b[i__2].i = z__1.i;
1269 /* Compute rhs (C, F) */
1271 zgemm_("N", "N", m, n, m, &c_b1, &a[a_offset], lda, &r__[r_offset], ldr, &
1272 c_b3, &c__[c_offset], ldc);
1273 z__1.r = -1., z__1.i = 0.;
1274 zgemm_("N", "N", m, n, n, &z__1, &l[l_offset], ldl, &b[b_offset], ldb, &
1275 c_b1, &c__[c_offset], ldc);
1276 zgemm_("N", "N", m, n, m, &c_b1, &d__[d_offset], ldd, &r__[r_offset], ldr,
1277 &c_b3, &f[f_offset], ldf);
1278 z__1.r = -1., z__1.i = 0.;
1279 zgemm_("N", "N", m, n, n, &z__1, &l[l_offset], ldl, &e[e_offset], lde, &
1280 c_b1, &f[f_offset], ldf);