/*
   Copyright (C) 2011  Statoil ASA, Norway. 
    
   The file 'matrix.h' is part of ERT - Ensemble based Reservoir Tool. 
    
   ERT 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 3 of the License, or 
   (at your option) any later version. 
    
   ERT 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 at <http://www.gnu.org/licenses/gpl.html> 
   for more details. 
*/

#ifndef __MATRIX_H__
#define __MATRIX_H__
#include <stdlib.h>
#include <stdio.h>
#include <stdbool.h>

#include <ert/util/rng.h>
#include <ert/util/type_macros.h>
#ifdef HAVE_THREAD_POOL
#include <ert/util/thread_pool.h>
#endif

#ifdef __cplusplus 
extern "C" {
#endif

typedef struct matrix_struct matrix_type;

  void          matrix_fscanf_data( matrix_type * matrix , bool row_major_order , FILE * stream );
  void          matrix_fprintf( const matrix_type * matrix , const char * fmt , FILE * stream );
  void          matrix_pretty_fprint(const matrix_type * matrix , const char * name , const char * fmt , FILE * stream);
  matrix_type * matrix_alloc(int rows, int columns);
  matrix_type * matrix_safe_alloc(int rows, int columns);
  bool          matrix_resize(matrix_type * matrix , int rows , int columns , bool copy_content);
  bool          matrix_safe_resize(matrix_type * matrix , int rows , int columns , bool copy_content);
  matrix_type * matrix_alloc_copy(const matrix_type * src);
  matrix_type * matrix_safe_alloc_copy(const matrix_type * src);
  matrix_type * matrix_realloc_copy(matrix_type * T , const matrix_type * src);
  
  matrix_type * matrix_alloc_shared(const matrix_type * src , int row , int column , int rows , int columns);
  void          matrix_free(matrix_type * matrix);
  void          matrix_safe_free( matrix_type * matrix );
  void          matrix_pretty_print(const matrix_type * matrix , const char * name , const char * fmt);
  void          matrix_set(matrix_type * matrix, double value);
  void          matrix_set_name( matrix_type * matrix , const char * name);
  void          matrix_scale(matrix_type * matrix, double value);
  void          matrix_shift(matrix_type * matrix, double value);
  
  void          matrix_assign(matrix_type * A , const matrix_type * B);
  void          matrix_inplace_add(matrix_type * A , const matrix_type * B);
  void          matrix_inplace_sub(matrix_type * A , const matrix_type * B);
  void          matrix_inplace_mul(matrix_type * A , const matrix_type * B);
  void          matrix_inplace_div(matrix_type * A , const matrix_type * B);
  void          matrix_sub(matrix_type * A , const matrix_type * B , const matrix_type * C);
  void          matrix_mul( matrix_type * A , const matrix_type * B , const matrix_type * C);
  void          matrix_transpose(const matrix_type * A , matrix_type * T);
  
  void          matrix_iset_safe(matrix_type * matrix , int i , int j, double value);
  void          matrix_iset(matrix_type * matrix , int i , int j, double value);
  double        matrix_iget(const matrix_type * matrix , int i , int j);
  double        matrix_iget_safe(const matrix_type * matrix , int i , int j);
  void          matrix_iadd(matrix_type * matrix , int i , int j , double value);
  void          matrix_isub(matrix_type * matrix , int i , int j , double value);
  void          matrix_imul(matrix_type * matrix , int i , int j , double value);
  
  
  void          matrix_inplace_matmul(matrix_type * A, const matrix_type * B);
  void          matrix_inplace_matmul_mt1(matrix_type * A, const matrix_type * B , int num_threads);
#ifdef HAVE_THREAD_POOL
  void          matrix_inplace_matmul_mt2(matrix_type * A, const matrix_type * B , thread_pool_type * thread_pool);
#endif
  
  void          matrix_shift_column(matrix_type * matrix , int column, double shift);
  void          matrix_shift_row(matrix_type * matrix , int row , double shift);
  double        matrix_get_column_sum(const matrix_type * matrix , int column);
  double        matrix_get_row_sum(const matrix_type * matrix , int column);
  double        matrix_get_column_sum2(const matrix_type * matrix , int column);
  double        matrix_get_row_abssum(const matrix_type * matrix , int row);
  double        matrix_get_column_abssum(const matrix_type * matrix , int column);
  double        matrix_get_row_sum2(const matrix_type * matrix , int column);
  void          matrix_subtract_row_mean(matrix_type * matrix);
  void          matrix_subtract_and_store_row_mean(matrix_type * matrix, matrix_type * row_mean);
  void          matrix_scale_column(matrix_type * matrix , int column  , double scale_factor);
  void          matrix_scale_row(matrix_type * matrix , int row  , double scale_factor);
  void          matrix_set_const_column(matrix_type * matrix , const double value , int column);
  void          matrix_copy_column(matrix_type * target_matrix, const matrix_type * src_matrix , int src_column, int target_column);
  void          matrix_set_const_row(matrix_type * matrix , const double value , int row);
  
  double      * matrix_get_data(const matrix_type * matrix);
  double        matrix_orthonormality( const matrix_type * matrix );
  
  matrix_type * matrix_alloc_steal_data(int rows , int columns , double * data , int data_size);
  void          matrix_set_column(matrix_type * matrix , const double * data , int column);
  void          matrix_set_many_on_column(matrix_type * matrix , int row_offset , int elements , const double * data , int column);
  void          matrix_ensure_rows(matrix_type * matrix, int rows, bool copy_content);
  void          matrix_shrink_header(matrix_type * matrix , int rows , int columns);
  void          matrix_full_size( matrix_type * matrix );
  int           matrix_get_rows(const matrix_type * matrix);
  int           matrix_get_columns(const matrix_type * matrix);
  int           matrix_get_row_stride(const matrix_type * matrix);
  int           matrix_get_column_stride(const matrix_type * matrix);
  void          matrix_get_dims(const matrix_type * matrix ,  int * rows , int * columns , int * row_stride , int * column_stride);
  bool          matrix_is_quadratic(const matrix_type * matrix);
  bool          matrix_equal( const matrix_type * m1 , const matrix_type * m2);
  
  void          matrix_diag_set_scalar(matrix_type * matrix , double value);
  void          matrix_diag_set(matrix_type * matrix , const double * diag);
  void          matrix_random_init(matrix_type * matrix , rng_type * rng);
  void          matrix_matlab_dump(const matrix_type * matrix, const char * filename);
  
  void          matrix_imul_col( matrix_type * matrix , int column , double factor);
  double        matrix_column_column_dot_product(const matrix_type * m1 , int col1 , const matrix_type * m2 , int col2);
  double        matrix_row_column_dot_product(const matrix_type * m1 , int row1 , const matrix_type * m2 , int col2);
  matrix_type * matrix_alloc_view(double * data , int rows , int columns);
  matrix_type * matrix_alloc_transpose( const matrix_type * A);
  void          matrix_copy_row(matrix_type * target_matrix, const matrix_type * src_matrix , int target_row, int src_row);
  void          matrix_copy_block( matrix_type * target_matrix , int target_row , int target_column , int rows , int columns,
                                   const matrix_type * src_matrix , int src_row , int src_column);
  
  void          matrix_scalar_set( matrix_type * matrix , double value);
  void          matrix_inplace_diag_sqrt(matrix_type *Cd);
  void          matrix_create_identiy(int n,matrix_type *Id);
  double        matrix_trace(const matrix_type *matrix);
  double        matrix_diag_std(const matrix_type * Sk,double mean);
  double        matrix_det3( const matrix_type * A);
  double        matrix_det4( const matrix_type * A);
  
  #ifdef HAVE_ISFINITE
  bool          matrix_is_finite(const matrix_type * matrix);
  void          matrix_assert_finite( const matrix_type * matrix );
  #endif

  UTIL_SAFE_CAST_HEADER( matrix );

#ifdef __cplusplus 
}
#endif
#endif
