/* a bug in zgelsd_ */

/*

setenv LD_LIBRARY_PATH /Components/MKL/11.0.2/Linux-x86-64/vendor

gcc -Wall -g -O0 -o zgelsd_bug zgelsd_bug.c -fno-strict-aliasing -L $LD_LIBRARY_PATH -lmkl_intel_ilp64 -lmkl_intel_thread -lmkl_core -liomp5 -lm

*/


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

/* from f2c.h */
typedef long int integer;
typedef double doublereal;
typedef struct { doublereal r, i; } doublecomplex;

/* zgelsd prototype */
int zgelsd_(integer *m, integer *n, integer *nrhs,
        doublecomplex *a, integer *lda, doublecomplex *b, integer *ldb,
        doublereal *s, doublereal *rcond, integer *rank, doublecomplex *work,
        integer *lwork, doublereal *rwork, integer *iwork, integer *info);

#define DIM 36

#define INPUTFNAME    "mat_rhs_vals.h"

void read_input_data(doublecomplex *A, doublecomplex *b, integer n)
{
   FILE *file;
   int i, j, k, l = 0, index;
   double tmp;
   
   file = fopen(INPUTFNAME, "r");
   if (file == NULL) {
      printf("The file %s cannot be found.\n", INPUTFNAME);
      exit(1);
   }

   printf("Reading the input matrix... \n");
   for (i = 0; i < n; i++) {
      for (j = 0; j < n; j++) {
         index = j*n+i; /* A is transposed */
         k = fscanf(file, "(%lg %lg) ", &A[index].r, &A[index].i);
         l++;
         if (k != 2 || k == EOF) {
            printf("\nWrong input data. Check the file %s.\n", INPUTFNAME);
            printf("Element # = %d, A[%d] = (%lg, %lg), n = %ld\n",
                   l, index, A[index].r, A[index].i, n);
            exit(1);
        }
      }
   }
   
   printf("Reading the input RHS... \n");
   for (i = 0; i < n; i++) {
      k = fscanf(file, "%lg ", &b[i].r);
      b[i].i = 0.0;
      if (k != 1 || k == EOF) {
         printf("\nWrong input data. Check the file %s.\n", INPUTFNAME);
         printf("Element # = %d, rhs[%d] = (%lg %lg), n = %ld\n", 
                i, i, b[i].r, b[i].i, n);
         exit(1);
      }
   }

   k = fscanf(file, "%lg ", &tmp);
   if (k == EOF)
      printf("Done\n");
   else 
      printf("The end of file is not reached...\n");
   fclose(file);
}


int main()
{
   integer m = DIM, n = DIM, nrhs = 1;
   integer mn_max, mn_min, rank, info;
   integer lwork, liwork, lrwork;
   integer smlsize = 25, nlvl, tmp1, tmp2;
   integer *iwork;
   doublecomplex A[DIM*DIM];
   doublecomplex B[DIM];
   doublecomplex *work, work1;
   doublereal S[DIM], *rwork;
   doublereal rcond;

   read_input_data(A, B, n);
   mn_max = (m > n) ? m : n;
   mn_min = (m > n) ? n : m;

   nlvl = (int) (log(((double)mn_min)/(smlsize+1.0))/log(2.0))+1;
   nlvl = (nlvl < 0) ? 0 : nlvl;          
   liwork = 3*mn_min*nlvl + 11*mn_min;
   iwork = (integer *) malloc(liwork*sizeof(integer));

   tmp1 = (smlsize+1)*(smlsize+1);
   tmp2 = n*(1+nrhs)+2*nrhs;
   lrwork = tmp1 > tmp2 ? tmp1 : tmp2;
   lrwork += 10*mn_min + 2*mn_min*smlsize + 8*mn_min*nlvl + 3*smlsize*nrhs;
   rwork = (doublereal *) malloc(lrwork*sizeof(doublereal));

   lwork = -1;
   zgelsd_(&m, &n, &nrhs, A, &m, B, &mn_max, 
           S, &rcond, &rank, &work1, &lwork, rwork, iwork, &info);
   lwork = (int) work1.r;
   work = (doublecomplex *) malloc(lwork*sizeof(doublecomplex));

   rcond = 100*2.2204460492503130808e-16;
   zgelsd_(&m, &n, &nrhs, A, &m, B, &mn_max, 
           S, &rcond, &rank, work, &lwork, rwork, iwork, &info);

   printf("info = %ld\n", info);

   return 0;
}