// ----------------------------------------------------------------------------
//
// Copyright (C) Columbia University, 1994. All Rights Reserved.
// Sameer A. Nene, Shree K. Nayar, Hiroshi Murase
//
// See file LICENSE for details of software license agreement.
//
// ----------------------------------------------------------------------------
//
// gso.C
//
// Author:               Sameer Nene
// Date:                 06/09/94
// Version:              1.0.0
// Modification History:
//
// Classes:
//
// Notes:
//   This command is used to compute a space which is orthogonal to the
//   supplied eigen spaces. See man pages for more details
//   
// ----------------------------------------------------------------------------

#include <iostream.h>
#include <math.h>
#include <stdlib.h>
#include "errorscope.h"
#include "vector.h"
#include "stringary.h"
#include "parser.h"
#include "ioutil.h"

#define GSO__ACCURACY .000001

main(char argc, char **argv)
{
  int i, j, k, l, m, n, o, p;
  double e;
  CommandParser cp(argc, argv, "e:V");
  ErrorScope::Error error;

  if((error = cp.error()) != ErrorScope::OK) {
    cerr << error << endl << "Exiting" << endl;
    exit(1);
  }

  if(cp.getVerboseFlag() != 0)
    cerr << "Copyright \xA9 Columbia University, 1994. All Rights Reserved." << endl << "Sameer Nene, Shree Nayar, Hiroshi Murase" << endl << endl << "Orthogonalizing..." << endl;

  StringAry &sa = cp.getObjectList();

  VectAry **va = new VectAry*[j = sa.getSize()];

  for(i = 0, k = 0; i < j; ++i) {
    va[i] = new VectAry;

    if((error = IOUtil::get(va[i], sa[i])) != ErrorScope::OK) {
      cerr << error << ": " << sa[i] << endl << "Exiting" << endl;
      exit(1);
    }
    
    k += va[i] -> getSize();
  }

//
// Gram-Schmidt Orthogonalization
//
// MGSCI Algorithm taken from:
// "Iterative Algorithms for Gram-Schmidt Orthogonalization" W. Hoffman
// Computing 41 1989, pp 335-348
//
  
  VectAry a(k), q(k);

  for(i = 0, n = 0; i < j; ++i) {
    for(l = 0, m = va[i] -> getSize(); l < m; ++l, ++n) {
      a[n] = (*(va[i]))[l];
      for(o = 0, q[n].length(p = a[n].length()); o < p; ++o)
	q[n].element(o, 0);
    }
    delete va[i];
  }

  delete va;

  for(i = 0; i < k; ++i) {
    Vector &t = a[i];
    e = 0.;
    do {
      for(j = 0; j < i; ++j)
	t -= q[j] * (q[j] * t);
      for(m = 0; m < i; ++m)
	if(fabs(e = t * q[m]) > GSO__ACCURACY)
	  break;
    }
    while(fabs(e) > GSO__ACCURACY);
    q[i] = t / t();
  }

  if((error = IOUtil::put(q, cp.getEigenVectorFileName())) != ErrorScope::OK) {
    cerr << error << ": " << cp.getEigenVectorFileName() << endl << "Exiting" << endl;
    exit(1);
  }

  if(cp.getVerboseFlag() != 0)
    cerr << "Done..." << endl;

  return 0;
}
