#include <stdlib.h>
#include <stdio.h>
#include "dgesvd_143t.h"


/* Parameters */
#define M MM
#define N NN
#define LDA M
#define LDU M
#define LDVT 1

/* Main program */
int main() {
    /* Locals */
    int m = M, n = N, lda = LDA, ldu = LDU, ldvt = LDVT, info, lwork;
    double wkopt;
    double* work;
    /* Local arrays */
/* iwork dimension should be at least 8*min(m,n) */
    int iwork[8 * N];
    double s[N], u[LDU * M], vt[LDVT * N];

    lwork = -1;
 
    dgesdd("N", &m, &n, a, &lda, s, u, &ldu, vt, &ldvt, &wkopt,
        &lwork, iwork, &info);
    lwork = (int)wkopt;
    work = (double*)malloc(lwork * sizeof(double));
    /* Compute SVD */
    dgesdd("N", &m, &n, a, &lda, s, u, &ldu, vt, &ldvt, work,
        &lwork, iwork, &info);
    /* Check for convergence */
    if (info > 0) {
        printf("The algorithm computing SVD failed to converge.\n");
        exit(1);
    }
    /* Print singular values */
    print_matrix("Singular values", 1, n, s, 1);

    /* Free workspace */
    free((void*)work);
    exit(0);
} /* End of DGESDD Example */

/* Auxiliary routine: printing a matrix */
void print_matrix(char* desc, int m, int n, double* a, int lda) {
    int i, j;
    printf("\n %s\n", desc);
    for (i = 0; i < m; i++) {
        for (j = 0; j < n; j++) 
            printf(" %le", a[i + j ]);
        printf("\n");
    }
}