/* nag_dbdsdc (f08mdc) Example Program. * * Copyright 2011 Numerical Algorithms Group. * * Mark 23, 2011. */ #include #include #include #include #include #include #include #include int main(void) { /* Scalars */ double alpha, beta, eps, norm; Integer abi, i, j, k1, k2, leniq, lenq, mlvl, n, pdab, pdb, pdu, pdvt; Integer exit_status = 0, smlsiz = 25; /* Arrays */ double *ab = 0, *b = 0, *d = 0, *e = 0, *q = 0, *u = 0, *vt = 0; Integer *iq = 0; char nag_enum_arg[40]; /* Nag Types */ NagError fail; Nag_OrderType order; Nag_UploType uplo; Nag_ComputeSingularVecsType compq; #ifdef NAG_COLUMN_MAJOR #define B(I, J) b[(J - 1) * pdb + I - 1] #define U(I, J) u[(J - 1) * pdu + I - 1] order = Nag_ColMajor; #else #define B(I, J) b[(I - 1) * pdb + J - 1] #define U(I, J) u[(I - 1) * pdu + J - 1] order = Nag_RowMajor; #endif INIT_FAIL(fail); printf("nag_dbdsdc (f08mdc) Example Program Results\n\n"); /* Skip heading in data file */ scanf("%*[^\n]"); scanf("%ld%*[^\n]", &n); if (n < 0) { printf("Invalid n\n"); exit_status = 1; goto END;; } scanf(" %s%*[^\n]", nag_enum_arg); /* nag_enum_name_to_value(x04nac). * Converts NAG enum member name to value */ uplo = (Nag_UploType) nag_enum_name_to_value(nag_enum_arg); scanf(" %s%*[^\n]", nag_enum_arg); /* Starting index for main diagonal in banded storage format = abi. */ if ((order==Nag_ColMajor && uplo==Nag_Lower) || (order==Nag_RowMajor && uplo==Nag_Upper)) { abi = 0; } else { abi = 1; } compq = (Nag_ComputeSingularVecsType) nag_enum_name_to_value(nag_enum_arg); /* size of u, vt, q and iq depends on value of compq input */ if (compq==Nag_SingularVecs) { pdu = n; pdvt = n; } else { pdu = 1; pdvt = 1; } if (compq==Nag_PackedSingularVecs) { mlvl = (Integer) (log(n/(smlsiz+1.0))/log(2.0)) + 1; if (mlvl < 1) mlvl = 1; lenq = MAX(n*n+5*n,n*(3 + 2*smlsiz + 8*mlvl)); leniq = n*3*mlvl; } else { lenq = 1; leniq = 1; } pdb = n; pdab = 2; /* Allocate memory */ if (!(b = NAG_ALLOC(n*n, double)) || !(ab = NAG_ALLOC(pdab*n, double)) || !(d = NAG_ALLOC(n, double)) || !(e = NAG_ALLOC(n-1, double)) || !(q = NAG_ALLOC(lenq, double)) || !(u = NAG_ALLOC(pdu*pdu, double)) || !(vt = NAG_ALLOC(pdvt*pdvt, double)) || !(iq = NAG_ALLOC(leniq, Integer))) { printf("Allocation failure\n"); exit_status = -1; goto END; } /* Read the bidiagonal matrix B from data file, * first the diagonal elements, and then the off-diagonal elements. */ for (i = 0; i < n; ++i) scanf("%lf", &d[i]); scanf("%*[^\n]"); for (i = 0; i < n - 1; ++i) scanf("%lf", &e[i]); scanf("%*[^\n]"); /* Store diagonal arrays in banded format in ab for printing */ for(i=0; i pow(eps,0.8)) { printf("Norm of B-(U*S*V^T) is much greater than 0.\nSchur " "factorization has failed.\n norm = %11.4e\n", norm); } } END: if (ab) NAG_FREE(ab); if (b) NAG_FREE(b); if (d) NAG_FREE(d); if (e) NAG_FREE(e); if (q) NAG_FREE(q); if (u) NAG_FREE(u); if (vt) NAG_FREE(vt); if (iq) NAG_FREE(iq); return exit_status; }