/***************************************************************************

  Parallel Linear Algebra Package Release R2.0.2

  6 Feb 2000
  
  Copyright (c) 1997,1998,1999,2000 Robert van de Geijn and 
  The University of Texas at Austin.
  See the file README for details on the gnu license

  This program is free software; you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation; either version 1, or (at your option)
  any later version.

  This program is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with this program; if not, write to the Free Software
  Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.

  Written under the direction of: 
 	Robert van de Geijn, Department of Computer Sciences,
	University of Texas at Austin  78712.    
	rvdg@cs.utexas.edu
 
  Many people contributed to the Parallel Linear Algebra Package (PLAPACK)
  including:	
  
  Gregory A. Baker, Center for Space Research,
  University of Texas at Austin  78722      baker@csr.utexas.edu

  Philip A. Alpatov, Department of Physics,
  University of Texas at Austin 78712.     philip@physics.utexas.edu

  James Overfelt, Department of Mathematics,
  University of Texas at Austin  78712.    overfelt@math.utexas.edu
  
  John Andrew Gunnels, Department of Computer Sciences,
  University of Texas at Austin 78712.     gunnels@cs.utexas.edu
  
  Greg Morrow, Department of Physics, 
  University of Texas at Austin 78712.     morrow@physics.utexas.edu

  Wesley Reiley, Department of Computer Sciences
  University of Texas at Austin 78712.     wesley@cs.utexas.edu

  Dr. Robert Van de Geijn, Department of Computer Sciences,
  University of Texas at Austin 78712.     rvdg@cs.utexas.edu

  This release (Release 2.0.2) was written primarily by
  Robert van de Geijn
 			
***********i****************************************************************/


#include "PLA.h"

static int old_size_malloced;

static PLA_Obj A_cpy = NULL;

int PLA_Chol_enter( int uplo, PLA_Obj A )
{
  int 
    value = PLA_SUCCESS,
    length_A, width_A, 
    objtype;
  
  MPI_Datatype
    datatype;

  char 
    routine_name[ 35 ] = "PLA_Chol";

  PLA_Routine_stack_push( routine_name );

  PLA_Routine_stack_push( "PLA_Chol_enter" );

  old_size_malloced = PLA_Total_size_malloced( );

  if ( PLA_CHECK_PARAMETERS ){
    /* Check if uplo parameter is valid */
    if ( !PLA_Valid_uplo_parameter( uplo ) ){
      PLA_Warning( "Invalid parameter uplo" );
      value--;
    }

    if ( uplo != PLA_LOWER_TRIANGULAR ){
      PLA_Abort( "Only lower triangular factorization supported",
		  __LINE__, __FILE__ );
      value--;
    }

    /* Check if A is valid matrix */

    if ( A == NULL || !PLA_Valid_object( A ) ) {
      PLA_Warning( "Invalid object A" );
      value--;
    }

    PLA_Obj_objtype( A, &objtype );
    if ( objtype != PLA_MATRIX ){
      PLA_Warning( "Invalid objtype for A" );
      value--;
    }      

    /* Check if matrix dimensions match */
    PLA_Obj_global_length( A, &length_A );
    PLA_Obj_global_width(  A, &width_A );

    if ( length_A != width_A ){
      PLA_Warning( "A is not square" );
      value--;
    }      

    /* Check if datatype is currently supported */
    PLA_Obj_datatype( A, &datatype );
    if ( datatype != MPI_FLOAT && datatype != MPI_DOUBLE 
	 && datatype != MPI_COMPLEX && datatype != MPI_DOUBLE_COMPLEX ){
      PLA_Warning( "datatype not yet supported" );
      value--;
    }
  }

  if ( PLA_CHECK_AGAINST_SEQUENTIAL ){
    PLA_Mscalar_create_conf_to( 
              A, PLA_ALL_ROWS, PLA_ALL_COLS, &A_cpy );

    PLA_Copy( A, A_cpy );
  }

  PLA_Routine_stack_pop( routine_name );
				 
  return value;
}

#define max(x,y) ( (x) > (y) ? (x) : (y) )

int PLA_Chol_exit( int uplo, PLA_Obj A )

{
  int value = PLA_SUCCESS,
      size_malloced;
  double 
    max_A, diff,
    PLA_Local_abs_max(), PLA_Local_abs_diff();
  char 
    routine_name[ 35 ];


  PLA_Routine_stack_push( "PLA_Chol_exit" );

  if ( PLA_CHECK_AGAINST_SEQUENTIAL ){
    PLA_Obj A_temp = NULL;
    MPI_Datatype
      datatype;
    double tol;

    PLA_Obj_datatype( A, &datatype );
    if ( datatype == MPI_DOUBLE || datatype == MPI_DOUBLE_COMPLEX )
      tol = 0.000001;
    else
      tol = 0.0001;

    PLA_Mscalar_create_conf_to( A, PLA_ALL_ROWS, PLA_ALL_COLS, &A_temp );
    PLA_Copy( A, A_temp );

    max_A = PLA_Local_abs_max( A_cpy );

    PLA_Local_chol( uplo, A_cpy );

    diff = PLA_Local_abs_diff( A_temp, A_cpy );

    if ( diff > tol * max_A ){
      PLA_Warning( "PLA_Chol: large relative error encountered" );
      value--;
    }      

    PLA_Obj_free( &A_cpy );
    PLA_Obj_free( &A_temp );
  }

  size_malloced = PLA_Total_size_malloced( );

  if ( size_malloced != old_size_malloced )
    PLA_Warning( "PLA_Chol: memory discrepency" );
  
  PLA_Routine_stack_pop( routine_name );

  PLA_Routine_stack_pop( routine_name );

  return value;
}