Macaulay2 Engine
Loading...
Searching...
No Matches

◆ cond_number_via_svd()

bool cond_number_via_svd ( int size,
complex * A,
double & cond )

Definition at line 1269 of file NAG.cpp.

1270{
1271 bool ret = true;
1272 char doit = 'A'; // other options are 'S' and 'O' for singular vectors only
1273 int rows = size;
1274 int cols = size;
1275 int info;
1276 int min = (rows <= cols) ? rows : cols;
1277
1278 if (min == 0)
1279 {
1280 ERROR("expected a matrix with positive dimensions");
1281 return false;
1282 }
1283
1284 int max = (rows >= cols) ? rows : cols;
1285 int wsize = 4 * min + 2 * max;
1286 double* workspace = newarray_atomic(double, 2 * wsize);
1287 double* rwork = newarray_atomic(double, 5 * max);
1288
1289 double* copyA = (double*)A;
1290 double* u = newarray_atomic(double, 2 * rows * rows);
1291 double* vt = newarray_atomic(double, 2 * cols * cols);
1292 double* sigma = newarray_atomic(double, 2 * min);
1293
1294 zgesvd_(&doit,
1295 &doit,
1296 &rows,
1297 &cols,
1298 copyA,
1299 &rows,
1300 sigma,
1301 u,
1302 &rows,
1303 vt,
1304 &cols,
1305 workspace,
1306 &wsize,
1307 rwork,
1308 &info);
1309
1310 if (info < 0)
1311 {
1312 ERROR("argument passed to zgesvd had an illegal value");
1313 ret = false;
1314 }
1315 else if (info > 0)
1316 {
1317 ERROR("zgesvd did not converge");
1318 ret = false;
1319 }
1320 else
1321 {
1322 cond = sigma[0] / sigma[size - 1];
1323 // print_complex_matrix(size,copyA);
1324 // printf("(s_large=%lf, s_small=%lf)\n", sigma[0], sigma[size-1]);
1325 }
1326
1327 freemem(workspace);
1328 freemem(rwork);
1329 // freemem(copyA);
1330 freemem(u);
1331 freemem(vt);
1332 freemem(sigma);
1333
1334 return ret;
1335}
int zgesvd_(char *jobU, char *jobV, int *rows, int *cols, double *A, int *ldA, double *Sigma, double *U, int *ldU, double *VT, int *ldVT, double *w, int *lwork, double *rwork, int *info)
void freemem(void *s)
Definition m2-mem.cpp:103
const int ERROR
Definition m2-mem.cpp:55
const mpreal min(const mpreal &x, const mpreal &y)
Definition mpreal.h:2792
#define newarray_atomic(T, len)
Definition newdelete.hpp:91
#define max(a, b)
Definition polyroots.cpp:52

References ERROR, freemem(), max, newarray_atomic, and zgesvd_().

Referenced by PathTracker::track().