```/* nag_lapackeig_dbdsdc (f08mdc) Example Program.
*
* Copyright 2020 Numerical Algorithms Group.
*
* Mark 27.1, 2020.
*/

#include <math.h>
#include <nag.h>
#include <stdio.h>

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_lapackeig_dbdsdc (f08mdc) Example Program Results\n\n");

/* Skip heading in data file */
scanf("%*[^\n]");
scanf("%" NAG_IFMT "%*[^\n]", &n);
if (n < 0) {
printf("Invalid n\n");
exit_status = 1;
goto END;
;
}
scanf(" %39s%*[^\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(" %39s%*[^\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 < n; i++)
ab[2 * i + abi] = d[i];
for (i = 0; i < n - 1; i++)
ab[2 * i + abi + 1] = e[i];

/* k1 = lower bandwidth, k2 = upper bandwidth */
k1 = (uplo == Nag_Upper ? 0 : 1);
k2 = 1 - k1;

/* Print Bidiagonal Matrix B stored in ab.
* nag_file_print_matrix_real_band (x04cec).
* Print real packed banded matrix (easy-to-use)
*/
fflush(stdout);
nag_file_print_matrix_real_band(order, n, n, k1, k2, ab, 2, "Matrix B", 0,
&fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_file_print_matrix_real_band (x04cec).\n%s\n",
fail.message);
exit_status = 1;
goto END;
}

/* Calculate the singular values and left and right singular vectors of B
* using nag_lapackeig_dbdsdc (f08mdc).
*/
nag_lapackeig_dbdsdc(order, uplo, compq, n, d, e, u, pdu, vt, pdvt, q, iq,
&fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_lapackeig_dbdsdc (f08mdc).\n%s\n", fail.message);
exit_status = 1;
goto END;
}

printf("\nSingular values\n");
for (i = 0; i < n; ++i)
printf(" %8.4f%s", d[i], i % 8 == 7 ? "\n" : "");
printf("\n\n");

if (compq == Nag_SingularVecs) {
/* Reconstruct bidiagonal matrix from decomposition:
* first, U <- U*S, then Compute B = U*S*V^T.
* nag_blast_dgemm (f16yac).
*/
for (i = 1; i <= n; i++)
for (j = 1; j <= n; j++)
U(i, j) = U(i, j) * d[j - 1];
alpha = 1.0;
beta = 0.0;
nag_blast_dgemm(order, Nag_NoTrans, Nag_NoTrans, n, n, n, alpha, u, pdu, vt,
pdvt, beta, b, pdb, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_blast_dgemm (f16yac).\n%s\n", fail.message);
exit_status = 1;
goto END;
}

/* Subtract original bidiagonal matrix:
* this should give a matrix close to zero.
*/
for (i = 1; i <= n; i++)
B(i, i) -= ab[2 * i - 2 + abi];
if (uplo == Nag_Upper)
for (i = 1; i <= n - 1; i++)
B(i, i + 1) -= ab[2 * i - 1 + abi];
else
for (i = 1; i <= n - 1; i++)
B(i + 1, i) -= ab[2 * i - 1 + abi];

/* nag_blast_dge_norm (f16rac): Find norm of matrix B and print warning if
* it is too large.
*/
nag_blast_dge_norm(order, Nag_OneNorm, n, n, b, pdb, &norm, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_blast_dge_norm (f16rac).\n%s\n", fail.message);
exit_status = 1;
goto END;
}
/* Get the machine precision, using nag_machine_precision (x02ajc) */
eps = nag_machine_precision;
if (norm > pow(eps, 0.8)) {
printf("Norm of B-(U*S*V^T) is much greater than 0.\nSchur "
"factorization has failed.\n norm = %13.4e\n",
norm);
}
}
END:
NAG_FREE(ab);
NAG_FREE(b);
NAG_FREE(d);
NAG_FREE(e);
NAG_FREE(q);
NAG_FREE(u);
NAG_FREE(vt);
NAG_FREE(iq);

return exit_status;
}
```