/**
* @version 1.2, 04/03/2002 Quitar errores
*/
#ifndef QRSOLV_H
#define QRSOLV_H
//----------------------------------------------------------------------------
#include "../matrixdouble/squaredmatrix.h"
#include "../matrixdouble/vectordouble.h"
#include "rsolv.cpp"
#include <iostream>
using namespace std;


// (posada)

// Solves the set of n linear equations A x = b. a[1..n][1..n], c[1..n], and d[1..n] are
// input as the output of the routine qrdcmp and are not modified. b[1..n] is input as the
// right-hand side vector, and is overwritten with the solution vector on output.
// antes void qrsolv(double **a, int n, double c[], double d[], double b[])
void qrsolv(SquaredMatrix &a, int n,
		VectorDouble &c, VectorDouble &d, VectorDouble &b)
{
        // cout << "\n{ qrsolv";
	// { posada
        // Logica
	a.setLogic(1);
	b.setLogic(1);
	c.setLogic(1);
	d.setLogic(1);
	// No teclear, compatibilidad de versiones int n = a.getSize();
        // } posada

   // posada #include "rsolv.cpp" se necesita para correr rsolv
   //void rsolv(SquaredMatrix &a, int n, VectorDouble &d, VectorDouble &b);

   int i,j;
   double sum,tau;

   for (j=1;j<n;j++) {  // Form QT*b.
      for (sum=0.0,i=j;i<=n;i++)
         sum += a[i][j]*b[i];
      tau=sum/c[j];
      for (i=j;i<=n;i++)
         b[i] -= tau*a[i][j];
   }
   rsolv(a,n,d,b); // Solve R x = QT*b.
        // cout << "\n  } qrsolv";
}
/* (C) Copr. 1986-92 Numerical Recipes Software ;#. */
//----------------------------------------------------------------------------
#endif
// Fin------------------------------------------------------------------------