#include <stdio.h>
#include <mkl.h>

int ISEED[4] = {9, 10, 11, 12};

int main()
{
    int info = 0;
    const int n = 10;
    const int nrhs_base = 1;
    const int nrhs_max = 50;
    const double thresh = 0.001;

    double *D = (double*)malloc(n * sizeof(double));
    double *DL = (double*)malloc((n - 1) * sizeof(double));
    double *DU = (double*)malloc((n - 1) * sizeof(double));

    double *DF_BASE = (double*)malloc(n * sizeof(double));
    double *DLF_BASE = (double*)malloc((n - 1) * sizeof(double));
    double *DUF_BASE = (double*)malloc((n - 1) * sizeof(double));
    double *B_BASE = (double*)malloc(n * nrhs_base * sizeof(double));

    const int idist = 2;
    const int n1 = n - 1;
    dlarnv_(&idist, ISEED, &n, D);
    dlarnv_(&idist, ISEED, &n1, DL);
    dlarnv_(&idist, ISEED, &n1, DU);

    for (int j = 0; j < nrhs_base; j++)
         dlarnv_(&idist, ISEED, &n, &B_BASE[j * n]);

    const int ione = 1;
    dcopy(&n, D, &ione, DF_BASE, &ione);
    dcopy(&n1, DL, &ione, DLF_BASE, &ione);
    dcopy(&n1, DU, &ione, DUF_BASE, &ione);

    // func call
    dgtsv_(&n, &nrhs_base, DLF_BASE, DF_BASE, DUF_BASE, B_BASE, &n, &info);
    if (info != 0)
        printf("gtsv info error for n = %d, nrhs = %d, info = %d\n", n, nrhs_base, info);

    for (int nrhs = 1; nrhs <= nrhs_max; nrhs++) {
        double *DF = (double*)malloc(n * sizeof(double));
        double *DLF = (double*)malloc((n - 1) * sizeof(double));
        double *DUF = (double*)malloc((n - 1) * sizeof(double));
        double *BF = (double*)malloc(n * nrhs * sizeof(double));

        dcopy(&n, D, &ione, DF, &ione);
        dcopy(&n1, DL, &ione, DLF, &ione);
        dcopy(&n1, DU, &ione, DUF, &ione);
        for (int j = 0; j < nrhs; j++)
            dlarnv_(&idist, ISEED, &n, &BF[j * n]);

        dgtsv_(&n, &nrhs, DLF, DF, DUF, BF, &n, &info);
        if (info != 0)
            printf("gtsv info error for n = %d, nrhs = %d, info = %d\n", n, nrhs, info);

        for (int i = 0; i < n; i++)
            DF[i] -= DF_BASE[i];

        const int n2 = n - 2;
        for (int i = 0; i < n2; i++)
            DLF[i] -= DLF_BASE[i];

        for (int i = 0; i < n1; i++)
            DUF[i] -= DUF_BASE[i];

        double normD  = dnrm2(&n, DF, &ione);
        double normDL = dnrm2(&n2, DLF, &ione);
        double normDU = dnrm2(&n1, DUF, &ione);

        if (normD > thresh || normDL > thresh || normDU > thresh)
            printf("FAILED for n = %d, nrhs = %d: normD = %lf, normDL = %lf, normDU = %lf\n", n, nrhs, normD, normDL, normDU);
        else
            printf("PASSED for n = %d, nrhs = %d\n", n, nrhs);

        free(DF);
        free(DLF);
        free(DUF);
        free(BF);
    }
    free(D);
    free(DL);
    free(DU);

    free(DF_BASE);
    free(DLF_BASE);
    free(DUF_BASE);
    free(B_BASE);
}