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


// modificado Daniel Alba Cuellar (posada)
#include <math.h>
#define NRANSI
#include "nrutil.h"

//Constructs the QR decomposition of a[1..n][1..n]. The upper triangular matrix R is re-
//turned in the upper triangle of a, except for the diagonal elements of R which are returned in
//d[1..n]. The orthogonal matrix Q is represented as a product of n . 1 Householder matrices
//Q 1 : ::Q n . 1 , whereQ j =1 . uj
// uj=cj. The ith component of uj is zero for i = 1; :: : ; j . 1
//while the nonzero components are returned in a[i][j] for i = j; : : :; n. sing returns as
//true (1) if singularity is encountered during the decomposition, but the decomposition is still
//completed in this case; otherwise it returns false (0).

// antes int qrdcmp(double **a, int n, double *c, double *d, int *sing)
int qrdcmp(SquaredMatrix &a, int n, VectorDouble &c, VectorDouble &d)
{
	//cout << "\n{ qrdcmp";
	// { posada
        // Logica
	a.setLogic(1);
	c.setLogic(1);
	d.setLogic(1);
	// No teclear, compatibilidad de versiones int n = a.getSize();
        // } posada

   int i,j,k;
   // antes double scale=0.0,sigma,sum,tau;
   double scale,sigma,sum,tau;

   // antes *sing=0;
   int sing=0;
   for (k=1;k<n;k++) {
      scale=0.0;	// nuevo
      for (i=k;i<=n;i++)
         scale=FMAX(scale,fabs(a[i][k]));
      if (scale == 0.0) {             // Singular case.
         // antes *sing=1;
         sing=1;
         c[k]=d[k]=0.0;
      } else {                        // Form Q k and Q k . A.
         for (i=k;i<=n;i++)
            a[i][k] /= scale;
         for (sum=0.0,i=k;i<=n;i++)
            sum += pow(a[i][k],2.0);
         sigma=SIGN(sqrt(sum),a[k][k]);
         a[k][k] += sigma;
         c[k]=sigma*a[k][k];
         d[k] = -scale*sigma;
         for (j=k+1;j<=n;j++) {
            for (sum=0.0,i=k;i<=n;i++)
               sum += a[i][k]*a[i][j];
            tau=sum/c[k];
            for (i=k;i<=n;i++)
               a[i][j] -= tau*a[i][k];
         }
      }
   }
   d[n]=a[n][n];
   if (d[n] == 0.0)
	// antes *sing=1;
      sing=1;
	//cout << "\n  } qrdcmp";
   return sing;
}
#undef NRANSI
/* (C) Copr. 1986-92 Numerical Recipes Software ;#. */
//----------------------------------------------------------------------------
#endif
// Fin------------------------------------------------------------------------