/*!
 * @file svd.cpp
 * @brief svd object
 *
 * @date 2003.10.30
 *
 * Copyright 2003 Takashi Suehiro
 */
#include <stdio.h>
#include "svd.h"
/*!
 * @class SVD
 * @brief class for svd object
 *
 */

//! row = 1..m, m <= m_max,  column = 1..n
SVD::SVD(int m_max, int n) {
  mm_max = m_max;
  nn = n;
  anr = matrix(1, mm_max, 1, nn);
  unr = matrix(1, mm_max, 1, nn);
  vnr = matrix(1, nn, 1, nn);
  wnr = vector(1, nn);
  xnr = vector(1, nn);
  ynr = vector(1, mm_max);
  mm_cur = 0;
}

SVD::~SVD() {
  free_matrix(anr, 1, mm_max, 1, nn);
  free_matrix(unr, 1, mm_max, 1, nn);
  free_matrix(vnr, 1, nn, 1, nn);
  free_vector(wnr,1, nn);
  free_vector(xnr,1, nn);
  free_vector(ynr,1, mm_max);
}

//! row = 1..m m <= m_max
int SVD::set_m(int m) {
  if(m < 0) {
    mm_cur = 0;
    return -1;
  } else if(m > mm_max) {
    mm_cur = mm_max;
    return -1;
  } else {
    mm_cur = m;
    return 0;
  }
}

//! 0 <= idx < m_max, converted to (idx + 1)
//! j[i], i = 0..(n-1 or greater)
int SVD::put_row(int idx, double * jb) {
  int i;
  if(idx < 0) {
    return -1;
  } else if (idx >= mm_max) {
    return -1;
  } else {
    for(i=0;i<nn;i++) {
      anr[idx+1][i+1] = jb[i];
    }
    return 0;
  }
}
//! ys[i], i = 0..(m-1 or greater)
int SVD::put_ys(double * ys) {
  int i;
  for(i=0;i<mm_cur;i++) {
    ynr[i+1] = ys[i];
  }
  return 0;
}
int SVD::put_row_and_y(int idx, double y, double * jb) {
  int i;
  if(idx < 0) {
    return -1;
  } else if (idx >= mm_max) {
    return -1;
  } else {
    ynr[idx+1] = y;
    for(i=0;i<nn;i++) {
      anr[idx+1][i+1] = jb[i];
    }
    return 0;
  }
}
int SVD::do_svdcmp() {
  int i,j;
  for(i=1;i<=mm_cur;i++) {
    for (j=1;j<=nn;j++) {
      unr[i][j]=anr[i][j];
    }
  }
  if (mm_cur < nn) {
    for(i=mm_cur+1;i<=nn;i++) {
      for (j=1;j<=nn;j++) {
        unr[i][j]=0.0;
      }
    }
    svdcmp(unr,nn,nn,wnr,vnr);
  } else {
    svdcmp(unr,mm_cur,nn,wnr,vnr);
  }
  return 0;
}
int SVD::do_round_off(double lim) {
  double wmax,wmin;
  int j;
  wmax=0.0;
  for(j=1;j<=nn;j++) {
    if (wnr[j] > wmax) wmax=wnr[j];
  }
  wmin=wmax*lim;
  for(j=1;j<=nn;j++) {
    if (wnr[j] < wmin) wnr[j]=0.0;
  }
  return 0;
}
int SVD::do_svbksb() {
  int i;
  if (mm_cur < nn) {
    for(i=mm_cur+1;i<=nn;i++) {
      ynr[i]=0.0;
    }
    svbksb(unr,wnr,vnr,nn,nn,ynr,xnr);
  } else {
    svbksb(unr,wnr,vnr,mm_cur,nn,ynr,xnr);
  }
  return 0;
}
//! x[i], i = 0..(n-1 or greater)
int SVD::get_sv(double * w) {
  int i;
  for(i=0;i<nn;i++) {
    w[i] = wnr[i+1];
  }
  return 0;
}
int SVD::get_result(double * x) {
  int i;
  for(i=0;i<nn;i++) {
    x[i] = xnr[i+1];
  }
  return 0;
}
