/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
/*             ********   ***                                 SparseLib++    */
/*          *******  **  ***       ***      ***               v. 1.5c        */
/*           *****      ***     ******** ********                            */
/*            *****    ***     ******** ********              R. Pozo        */
/*       **  *******  ***   **   ***      ***                 K. Remington   */
/*        ********   ********                                 A. Lumsdaine   */
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
/*                                                                           */
/*                                                                           */
/*                     SparseLib++ : Sparse Matrix Library                   */
/*                                                                           */
/*               National Institute of Standards and Technology              */
/*                        University of Notre Dame                           */
/*              Authors: R. Pozo, K. Remington, A. Lumsdaine                 */
/*                                                                           */
/*                                 NOTICE                                    */
/*                                                                           */
/* Permission to use, copy, modify, and distribute this software and         */
/* its documentation for any purpose and without fee is hereby granted       */
/* provided that the above notice appear in all copies and supporting        */
/* documentation.                                                            */
/*                                                                           */
/* Neither the Institutions (National Institute of Standards and Technology, */
/* University of Notre Dame) nor the Authors make any representations about  */
/* the suitability of this software for any purpose.  This software is       */
/* provided ``as is'' without expressed or implied warranty.                 */
/*                                                                           */
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/

#include <stdlib.h>
#include "ilupre_double_GR.h" // modif GR : 
// 1) découle d'une classe virtuelle pour avoir une interface standard
// 2) introduction de deux fonctions d'initialisation qui servent pour les constructeurs
// initiaux et aussi pour les nouveaux constructeurs, ce qui évite la duplication
// de code
// 3) introduction de la class générique ILUPreconditioner_double, qui permet d'être
// au même niveau d'appel qu'avec les autres préconditionneurs

#include "qsort_double.h"

#include "spblas.h"
#include "Mat_creuse_CompCol.h"

// classe générique qui regroupe les différentes possibilités de 
// stockage matricielle 

// CONSTRUCTEUR

// cas d'une matrice abstraite
ILUPreconditioner_double::ILUPreconditioner_double(const Mat_abstraite &A) :
  preCond_mat(NULL)
   {
    switch (A.Type_matrice())
	 { case CREUSE_COMPRESSEE_COLONNE :
	    { preCond_mat = 
	        new CompCol_ILUPreconditioner_double(A);
		  break;
		 }
	   case CREUSE_COMPRESSEE_LIGNE :
	     // ici il y aura une modif a faire comme pour les matrices compressées colonne
	     // il faut passer la matrice dérivée et non la matrice mère
	    {cout << "\n problème dans ILUPreconditioner_double::ILUPreconditioner_double(.. "
	          << "voir le source .cpp, il y a du travail a faire";
	     Sortie(1);     
	     preCond_mat = 
	        new CompRow_ILUPreconditioner_double((CompRow_Mat_double&)A);
		  break;
		 }
	   default :
			cout << "\nErreur : valeur incorrecte du type de matrice !"
			     << " = " << A.Type_matrice() << " , ce tye n'est pas "
			     << " pris en compte dans le préconditionnement ILU pour cette classe\n";
			cout << "\n ILUPreconditioner_double:: "
			     << "ILUPreconditioner_double(... )\n";
			Sortie(1);
	  };
     };
       
// destructeur
ILUPreconditioner_double::~ILUPreconditioner_double()
  { if (preCond_mat != NULL) delete preCond_mat ;
   };

// ------- définition des classes spécifiques ------------- 

CompCol_ILUPreconditioner_double::
CompCol_ILUPreconditioner_double(const CompCol_Mat_double &A)
: 
  l_val_(0), l_colptr_(A.dim(1) + 1), l_rowind_(0), l_nz_(0), u_val_(0),
    u_colptr_(A.dim(1) + 1), u_rowind_(0), u_nz_(0)
      { Init_CompCol_Mat_double(A); 
        };
  
// cas d'une matrice abstraite
CompCol_ILUPreconditioner_double::
      CompCol_ILUPreconditioner_double(const Mat_abstraite &A) : 
  l_val_(0), l_colptr_(A.Nb_colonne() + 1), l_rowind_(0), l_nz_(0), u_val_(0),
    u_colptr_(A.Nb_colonne() + 1), u_rowind_(0), u_nz_(0)
 { 
   switch (A.Type_matrice())
	 { case CREUSE_COMPRESSEE_COLONNE :
	    {Init_Mat_creuse_CompCol(A);
		 break;
		 }
	   default :
			cout << "\nErreur : valeur incorrecte du type de matrice !"
			     << " = " << A.Type_matrice() << " , ce tye n'est pas "
			     << " pris en compte dans le préconditionnement ILU pour cette classe\n";
			cout << "\n CompCol_ILUPreconditioner_double:: "
			     << "CompCol_ILUPreconditioner_double(... )\n";
			Sortie(1);
	  };
   };

CompRow_ILUPreconditioner_double::
CompRow_ILUPreconditioner_double(const CompRow_Mat_double &A)
: l_val_(0), l_rowptr_(A.dim(1) + 1), l_colind_(0), l_nz_(0),
    u_val_(0), u_rowptr_(A.dim(1) + 1), u_colind_(0), u_nz_(0)
     { Init_CompRow_Mat_double( A);
       };

// cas d'une matrice abstraite
CompRow_ILUPreconditioner_double::CompRow_ILUPreconditioner_double
    (const Mat_abstraite &A) :
   l_val_(0), l_rowptr_(A.Nb_colonne() + 1), l_colind_(0), l_nz_(0),
    u_val_(0), u_rowptr_(A.Nb_colonne() + 1), u_colind_(0), u_nz_(0)
 { 
   switch (A.Type_matrice())
	 { case CREUSE_COMPRESSEE_LIGNE:
	    {Init_CompRow_Mat_double((CompRow_Mat_double&) A);
		 break;
		 }
	   default :
			cout << "\nErreur : valeur incorrecte du type de matrice !"
			     << " = " << A.Type_matrice() << " , ce tye n'est pas "
			     << " pris en compte dans le préconditionnement ILU pour cette classe\n";
			cout << "\n CompRow_ILUPreconditioner_double:: "
			     << "CompRow_ILUPreconditioner_double(... )\n";
			Sortie(1);
	  };
   };

VECTOR_double
CompCol_ILUPreconditioner_double::solve(const VECTOR_double &x) const 
{
  int M = x.size();
  VECTOR_double y(M);
  VECTOR_double work(M);

  int descra[9];

  descra[0] = 0;

  // lower unit
  descra[1] = 1;
  descra[2] = 1;

  F77NAME(dcscsm) (0, M, 1, 1, NULL, 1.0,
           descra, &l_val_(0), &l_rowind_(0), &l_colptr_(0),
           &x(0), M, 0.0, &y(1), M,
           &work(0), M);

  // upper diag
  descra[1] = 2;
  descra[2] = 0;

  F77NAME(dcscsm) (0, M, 1, 1, NULL, 1.0,
           descra, &u_val_(0), &u_rowind_(0), &u_colptr_(0),
           &y(0), M, 0.0, &y(1), M,
           &work(0), M);

  return y;
} 


VECTOR_double
CompCol_ILUPreconditioner_double::trans_solve(const VECTOR_double &x) const 
{
  int M = x.size();
  VECTOR_double y(M);
  VECTOR_double work(M);

  int descra[9];

  descra[0] = 0;

  // upper diag transpose
  descra[1] = 2;
  descra[2] = 0;

  F77NAME(dcscsm) (1, M, 1, 1, NULL, 1.0,
           descra, &u_val_(0), &u_rowind_(0), &u_colptr_(0),
           &x(0), M, 0.0, &y(1), M,
           &work(0), M);

  // lower unit transpose
  descra[1] = 1;
  descra[2] = 1;

  F77NAME(dcscsm) (1, M, 1, 1, NULL, 1.0,
           descra, &l_val_(0), &l_rowind_(0), &l_colptr_(0),
           &y(0), M, 0.0, &y(1), M,
           &work(0), M);

  return y;
}


VECTOR_double
CompRow_ILUPreconditioner_double::solve(const VECTOR_double &x) const 
{
  int M = x.size();
  VECTOR_double y(M);
  VECTOR_double work(M);

  int descra[9];

  descra[0] = 0;

  // lower unit
  descra[1] = 1;
  descra[2] = 1;

  F77NAME(dcsrsm) (0, M, 1, 1, NULL, 1.0,
           descra, &l_val_(0), &l_colind_(0), &l_rowptr_(0),
           &x(0), M, 0.0, &y(1), M,
           &work(0), M);

  // upper diag
  descra[1] = 2;
  descra[2] = 0;

  F77NAME(dcsrsm) (0, M, 1, 1, NULL, 1.0,
           descra, &u_val_(0), &u_colind_(0), &u_rowptr_(0),
           &y(0), M, 0.0, &y(1), M,
           &work(0), M);

  return y;
} 


VECTOR_double
CompRow_ILUPreconditioner_double::trans_solve(const VECTOR_double &x) const 
{
  int M = x.size();
  VECTOR_double y(M);
  VECTOR_double work(M);

  int descra[9];

  descra[0] = 0;

  // upper diag transpose
  descra[1] = 2;
  descra[2] = 0;

  F77NAME(dcsrsm) (1, M, 1, 1, NULL, 1.0,
           descra, &u_val_(0), &u_colind_(0), &u_rowptr_(0),
           &x(0), M, 0.0, &y(1), M,
           &work(0), M);

  // lower unit transpose
  descra[1] = 1;
  descra[2] = 1;

  F77NAME(dcsrsm) (1, M, 1, 1, NULL, 1.0,
           descra, &l_val_(0), &l_colind_(0), &l_rowptr_(0),
           &y(0), M, 0.0, &y(1), M,
           &work(0), M);

  return y;
} 
  
//====================== protégée ========================

  // fonction interne d'initialisation pour éviter la recopie
void CompCol_ILUPreconditioner_double::Init_CompCol_Mat_double
         (const CompCol_Mat_double &A)
 {
  int i, j, k, pn, qn, rn;
  double multiplier;

  // Copy
  dim_[0] = A.dim(0);
  dim_[1] = A.dim(1);


  // Get size of l and u
  for (i = 0; i < dim_[1]; i++)
    for (j = A.col_ptr(i); j < A.col_ptr(i+1); j++)
      if (A.row_ind(j) > i)
    l_nz_++;
      else
    u_nz_++;

  l_val_.newsize(l_nz_);
  u_val_.newsize(u_nz_);
  l_rowind_.newsize(l_nz_);
  u_rowind_.newsize(u_nz_);

  l_colptr_(0) = u_colptr_(0) = 0;

  // Split up A into l and u
  for (i = 0; i < dim_[1]; i++) {
    l_colptr_(i+1) = l_colptr_(i);
    u_colptr_(i+1) = u_colptr_(i);
    
    for (j = A.col_ptr(i); j < A.col_ptr(i+1); j++)
      if (A.row_ind(j) > i) {
    k = l_colptr_(i+1)++;
    l_val_(k) = A.val(j);
    l_rowind_(k) = A.row_ind(j);
      } else if (A.row_ind(j) <= i) {
    k = u_colptr_(i+1)++;
    u_val_(k) = A.val(j);
    u_rowind_(k) = A.row_ind(j);
      }
  }
  
  for (i = 0; i < dim_[1]; i++) {
    QSort(l_rowind_, l_val_, l_colptr_[i], l_colptr_[i+1] - l_colptr_[i]);
    QSort(u_rowind_, u_val_, u_colptr_[i], u_colptr_[i+1] - u_colptr_[i]);
  }

  // Factor matrix
  for (i = 0; i < dim_[0] - 1; i++) {
    multiplier = u_val_(u_colptr_(i+1)-1);
    
    for (j = l_colptr_[i]; j < l_colptr_[i+1]; j++)
      l_val_[j] /= multiplier;
    
    for (j = u_colptr_[i+1]; j < u_colptr_[i+2]-1; j++) {
      multiplier = u_val_[j];
      qn = j + 1;
      rn = l_colptr_[i+1];
      for (pn = l_colptr_[u_rowind_[j]]; 
       pn < l_colptr_[u_rowind_[j]+1] && l_rowind_[pn] <= i + 1; 
       pn++) {
    while (qn < u_colptr_[i+2] && u_rowind_[qn] < l_rowind_[pn])
      qn++;
    if (qn < u_colptr_[i+2] && l_rowind_[pn] == u_rowind_[qn])
      u_val_[qn] -= multiplier * l_val_[pn];
      }
      for ( ; pn < l_colptr_[u_rowind_[j]+1]; pn++) {
    while (rn < l_colptr_[i+2] && l_rowind_[rn] < l_rowind_[pn])
      rn++;
    if (rn < l_colptr_[i+2] && l_rowind_[pn] == l_rowind_[rn])
      l_val_[rn] -= multiplier * l_val_[pn];
      }
    }
  }
 };
  
void CompCol_ILUPreconditioner_double::Init_Mat_creuse_CompCol
 (const Mat_abstraite &B)
 {Mat_creuse_CompCol& A = (Mat_creuse_CompCol&) B;
  int i, j, k, pn, qn, rn;
  double multiplier;

  // Copy
  dim_[0] = A.dim(0);
  dim_[1] = A.dim(1);


  // Get size of l and u
  for (i = 0; i < dim_[1]; i++)
    for (j = A.col_ptr(i); j < A.col_ptr(i+1); j++)
      if (A.row_ind(j) > i)
    l_nz_++;
      else
    u_nz_++;

  l_val_.newsize(l_nz_);
  u_val_.newsize(u_nz_);
  l_rowind_.newsize(l_nz_);
  u_rowind_.newsize(u_nz_);

  l_colptr_(0) = u_colptr_(0) = 0;

  // Split up A into l and u
  for (i = 0; i < dim_[1]; i++) {
    l_colptr_(i+1) = l_colptr_(i);
    u_colptr_(i+1) = u_colptr_(i);
    
    for (j = A.col_ptr(i); j < A.col_ptr(i+1); j++)
      if (A.row_ind(j) > i) {
    k = l_colptr_(i+1)++;
    l_val_(k) = A.val(j);
    l_rowind_(k) = A.row_ind(j);
      } else if (A.row_ind(j) <= i) {
    k = u_colptr_(i+1)++;
    u_val_(k) = A.val(j);
    u_rowind_(k) = A.row_ind(j);
      }
  }
  
  for (i = 0; i < dim_[1]; i++) {
    QSort(l_rowind_, l_val_, l_colptr_[i], l_colptr_[i+1] - l_colptr_[i]);
    QSort(u_rowind_, u_val_, u_colptr_[i], u_colptr_[i+1] - u_colptr_[i]);
  }

  // Factor matrix
  for (i = 0; i < dim_[0] - 1; i++) {
    multiplier = u_val_(u_colptr_(i+1)-1);
    
    for (j = l_colptr_[i]; j < l_colptr_[i+1]; j++)
      l_val_[j] /= multiplier;
    
    for (j = u_colptr_[i+1]; j < u_colptr_[i+2]-1; j++) {
      multiplier = u_val_[j];
      qn = j + 1;
      rn = l_colptr_[i+1];
      for (pn = l_colptr_[u_rowind_[j]]; 
       pn < l_colptr_[u_rowind_[j]+1] && l_rowind_[pn] <= i + 1; 
       pn++) {
    while (qn < u_colptr_[i+2] && u_rowind_[qn] < l_rowind_[pn])
      qn++;
    if (qn < u_colptr_[i+2] && l_rowind_[pn] == u_rowind_[qn])
      u_val_[qn] -= multiplier * l_val_[pn];
      }
      for ( ; pn < l_colptr_[u_rowind_[j]+1]; pn++) {
    while (rn < l_colptr_[i+2] && l_rowind_[rn] < l_rowind_[pn])
      rn++;
    if (rn < l_colptr_[i+2] && l_rowind_[pn] == l_rowind_[rn])
      l_val_[rn] -= multiplier * l_val_[pn];
      }
    }
  }
 };
  
   
void CompRow_ILUPreconditioner_double::Init_CompRow_Mat_double
         (const CompRow_Mat_double &A)
{
  int i, j, k, pn, qn, rn;
  double multiplier;

  // Copy
  dim_[0] = A.dim(0);
  dim_[1] = A.dim(1);


  // Get size of l and u
  for (i = 0; i < dim_[1]; i++)
    for (j = A.row_ptr(i); j < A.row_ptr(i+1); j++)
      if (A.col_ind(j) < i)
    l_nz_++;
      else
    u_nz_++;

  l_val_.newsize(l_nz_);
  u_val_.newsize(u_nz_);
  l_colind_.newsize(l_nz_);
  u_colind_.newsize(u_nz_);

  l_rowptr_(0) = u_rowptr_(0) = 0;

  // Split up A into l and u
  for (i = 0; i < dim_[1]; i++) {
    l_rowptr_(i+1) = l_rowptr_(i);
    u_rowptr_(i+1) = u_rowptr_(i);
    
    for (j = A.row_ptr(i); j < A.row_ptr(i+1); j++)
      if (A.col_ind(j) < i) {
    k = l_rowptr_(i+1)++;
    l_val_(k) = A.val(j);
    l_colind_(k) = A.col_ind(j);
      } else if (A.col_ind(j) >= i) {
    k = u_rowptr_(i+1)++;
    u_val_(k) = A.val(j);
    u_colind_(k) = A.col_ind(j);
      }
  }

  for (i = 0; i < dim_[1]; i++) {
    QSort(l_colind_, l_val_, l_rowptr_[i], l_rowptr_[i+1] - l_rowptr_[i]);
    QSort(u_colind_, u_val_, u_rowptr_[i], u_rowptr_[i+1] - u_rowptr_[i]);
  }

  // Factor matrix
  for (i = 1; i < dim_[0]; i++) {
    for (j = l_rowptr_[i]; j < l_rowptr_[i+1]; j++) {
      pn = u_rowptr_[l_colind_[j]];
      multiplier = (l_val_[j] /= u_val_[pn]);

      qn = j + 1;
      rn = u_rowptr_[i];

      for (pn++; pn < u_rowptr_[l_colind_[j]+1] && u_colind_[pn] < i; pn++) {
    while (qn < l_rowptr_[i+1] && l_colind_[qn] < u_colind_[pn])
      qn++;
    if (qn < l_rowptr_[i+1] && u_colind_[pn] == l_colind_[qn])
      l_val_[qn] -= multiplier * u_val_[pn];
      }
      for ( ; pn < u_rowptr_[l_colind_[j]+1]; pn++) {
    while (rn < u_rowptr_[i+1] && u_colind_[rn] < u_colind_[pn])
      rn++;
    if (rn < u_rowptr_[i+1] && u_colind_[pn] == u_colind_[rn])
      u_val_[rn] -= multiplier * u_val_[pn];
      }
    }
  }
 };