00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013 #include "f2c.h"
00014 #include "blaswrap.h"
00015
00016
00017
00018 static complex c_b1 = {0.f,0.f};
00019 static complex c_b2 = {1.f,0.f};
00020
00021 int cget54_(integer *n, complex *a, integer *lda, complex *b,
00022 integer *ldb, complex *s, integer *lds, complex *t, integer *ldt,
00023 complex *u, integer *ldu, complex *v, integer *ldv, complex *work,
00024 real *result)
00025 {
00026
00027 integer a_dim1, a_offset, b_dim1, b_offset, s_dim1, s_offset, t_dim1,
00028 t_offset, u_dim1, u_offset, v_dim1, v_offset, i__1;
00029 real r__1, r__2;
00030 complex q__1;
00031
00032
00033 real dum[1], ulp, unfl;
00034 extern int cgemm_(char *, char *, integer *, integer *,
00035 integer *, complex *, complex *, integer *, complex *, integer *,
00036 complex *, complex *, integer *);
00037 real wnorm;
00038 extern doublereal clange_(char *, integer *, integer *, complex *,
00039 integer *, real *), slamch_(char *);
00040 extern int clacpy_(char *, integer *, integer *, complex
00041 *, integer *, complex *, integer *);
00042 real abnorm;
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140
00141 a_dim1 = *lda;
00142 a_offset = 1 + a_dim1;
00143 a -= a_offset;
00144 b_dim1 = *ldb;
00145 b_offset = 1 + b_dim1;
00146 b -= b_offset;
00147 s_dim1 = *lds;
00148 s_offset = 1 + s_dim1;
00149 s -= s_offset;
00150 t_dim1 = *ldt;
00151 t_offset = 1 + t_dim1;
00152 t -= t_offset;
00153 u_dim1 = *ldu;
00154 u_offset = 1 + u_dim1;
00155 u -= u_offset;
00156 v_dim1 = *ldv;
00157 v_offset = 1 + v_dim1;
00158 v -= v_offset;
00159 --work;
00160
00161
00162 *result = 0.f;
00163 if (*n <= 0) {
00164 return 0;
00165 }
00166
00167
00168
00169 unfl = slamch_("Safe minimum");
00170 ulp = slamch_("Epsilon") * slamch_("Base");
00171
00172
00173
00174 clacpy_("Full", n, n, &a[a_offset], lda, &work[1], n);
00175 clacpy_("Full", n, n, &b[b_offset], ldb, &work[*n * *n + 1], n)
00176 ;
00177
00178 i__1 = *n << 1;
00179 r__1 = clange_("1", n, &i__1, &work[1], n, dum);
00180 abnorm = dmax(r__1,unfl);
00181
00182
00183
00184 clacpy_(" ", n, n, &a[a_offset], lda, &work[1], n);
00185 cgemm_("N", "N", n, n, n, &c_b2, &u[u_offset], ldu, &s[s_offset], lds, &
00186 c_b1, &work[*n * *n + 1], n);
00187
00188 q__1.r = -1.f, q__1.i = -0.f;
00189 cgemm_("N", "C", n, n, n, &q__1, &work[*n * *n + 1], n, &v[v_offset], ldv,
00190 &c_b2, &work[1], n);
00191
00192
00193
00194 clacpy_(" ", n, n, &b[b_offset], ldb, &work[*n * *n + 1], n);
00195 cgemm_("N", "N", n, n, n, &c_b2, &u[u_offset], ldu, &t[t_offset], ldt, &
00196 c_b1, &work[(*n << 1) * *n + 1], n);
00197
00198 q__1.r = -1.f, q__1.i = -0.f;
00199 cgemm_("N", "C", n, n, n, &q__1, &work[(*n << 1) * *n + 1], n, &v[
00200 v_offset], ldv, &c_b2, &work[*n * *n + 1], n);
00201
00202
00203
00204 i__1 = *n << 1;
00205 wnorm = clange_("1", n, &i__1, &work[1], n, dum);
00206
00207 if (abnorm > wnorm) {
00208 *result = wnorm / abnorm / ((*n << 1) * ulp);
00209 } else {
00210 if (abnorm < 1.f) {
00211
00212 r__1 = wnorm, r__2 = (*n << 1) * abnorm;
00213 *result = dmin(r__1,r__2) / abnorm / ((*n << 1) * ulp);
00214 } else {
00215
00216 r__1 = wnorm / abnorm, r__2 = (real) (*n << 1);
00217 *result = dmin(r__1,r__2) / ((*n << 1) * ulp);
00218 }
00219 }
00220
00221 return 0;
00222
00223
00224
00225 }