OilfieldToolsNet/ResInsight
4
0
Fork 0
mirror of https://github.com/OPM/ResInsight.git synced 2025-02-25 18:55:39 -06:00
Code Issues Packages Projects Releases Wiki Activity

Upgraded ERT from github master commit d9955e849

Browse Source 
Startup on new project
This commit is contained in:
Jacob Støren
2015-04-14 15:47:22 +02:00
parent 02e47c6679
commit 09744af4b8
1119 changed files with 28508 additions and 12532 deletions
Download Patch File Download Diff File Expand all files Collapse all files

96
ThirdParty/Ert/devel/libanalysis/include/ert/analysis/analysis_module.h vendored
View File

@@ -1,19 +1,19 @@
/*
Copyright (C) 2011 Statoil ASA, Norway.
The file 'analysis_module.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.
Copyright (C) 2011 Statoil ASA, Norway.
The file 'analysis_module.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 __ANALYSIS_MODULE_H__
@@ -27,7 +27,7 @@ extern "C" {
#include <ert/util/bool_vector.h>
/*
/*
These are option flag values which are used by the core ert code to
query the module of it's needs and capabilities. For instance to to
determine whether the data should be scaled prior to analysis the
@@ -36,7 +36,7 @@ extern "C" {
if (analysis_module_get_option( module, ANALYSIS_SCALE_DATA))
obs_data_scale( obs_data , S , E , D , R , dObs );
It is the responsability of the module to set the various flags.
It is the responsability of the module to set the various flags.
*/
@@ -57,54 +57,54 @@ typedef enum {
{.value = ANALYSIS_ITERABLE , .name = "ANALYSIS_ITERABLE"}
#define EXTERNAL_MODULE_NAME "analysis_table"
#define EXTERNAL_MODULE_NAME "analysis_table"
#define EXTERNAL_MODULE_SYMBOL analysis_table
typedef enum {
LOAD_OK = 0,
DLOPEN_FAILURE = 1,
DLOPEN_FAILURE = 1,
LOAD_SYMBOL_TABLE_NOT_FOUND = 2
} analysis_module_load_status_enum;
typedef struct analysis_module_struct analysis_module_type;
analysis_module_type * analysis_module_alloc_internal__( rng_type * rng , const char * user_name , const char * symbol_table , bool verbose , analysis_module_load_status_enum * load_status);
analysis_module_type * analysis_module_alloc_internal( rng_type * rng , const char * user_name , const char * symbol_table );
analysis_module_type * analysis_module_alloc_external__(rng_type * rng , const char * user_name , const char * lib_name , bool verbose , analysis_module_load_status_enum * load_status);
analysis_module_type * analysis_module_alloc_external( rng_type * rng , const char * user_name , const char * libname );
void analysis_module_free( analysis_module_type * module );
void analysis_module_free__( void * arg);
void analysis_module_initX(analysis_module_type * module ,
matrix_type * X ,
matrix_type * A ,
matrix_type * S ,
matrix_type * R ,
matrix_type * dObs ,
matrix_type * E ,
void analysis_module_initX(analysis_module_type * module ,
matrix_type * X ,
matrix_type * A ,
matrix_type * S ,
matrix_type * R ,
matrix_type * dObs ,
matrix_type * E ,
matrix_type * D);
void analysis_module_updateA(analysis_module_type * module ,
matrix_type * A ,
matrix_type * S ,
matrix_type * R ,
matrix_type * dObs ,
matrix_type * E ,
void analysis_module_updateA(analysis_module_type * module ,
matrix_type * A ,
matrix_type * S ,
matrix_type * R ,
matrix_type * dObs ,
matrix_type * E ,
matrix_type * D );
void analysis_module_init_update( analysis_module_type * module ,
const bool_vector_type * ens_mask ,
const matrix_type * S ,
const matrix_type * R ,
const matrix_type * dObs ,
const matrix_type * E ,
void analysis_module_init_update( analysis_module_type * module ,
const bool_vector_type * ens_mask ,
const matrix_type * S ,
const matrix_type * R ,
const matrix_type * dObs ,
const matrix_type * E ,
const matrix_type * D );
const char * analysis_module_get_lib_name( const analysis_module_type * module);
bool analysis_module_internal( const analysis_module_type * module );
@@ -120,7 +120,7 @@ typedef enum {
bool analysis_module_get_bool( const analysis_module_type * module , const char * var);
void * analysis_module_get_ptr( const analysis_module_type * module , const char * var);
const char * analysis_module_flag_enum_iget( int index, int * value);
UTIL_IS_INSTANCE_HEADER( analysis_module );

2
ThirdParty/Ert/devel/libanalysis/modules/CMakeLists.txt vendored
View File

@@ -1,6 +1,8 @@
include_directories(${CMAKE_CURRENT_SOURCE_DIR})
set( RML_SOURCE_FILES
rml_enkf_config.c
rml_enkf_log.c
rml_enkf.c
rml_enkf_common.c )

518
ThirdParty/Ert/devel/libanalysis/modules/rml_enkf.c vendored
View File

@@ -1,19 +1,19 @@
/*
Copyright (C) 2011 Statoil ASA, Norway.
The file 'rml_enkf.c' 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.
Copyright (C) 2011 Statoil ASA, Norway.
The file 'rml_enkf.c' 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.
*/
@@ -35,6 +35,8 @@
#include <ert/analysis/std_enkf.h>
#include <rml_enkf_common.h>
#include <rml_enkf_config.h>
#include <rml_enkf_log.h>
typedef struct rml_enkf_data_struct rml_enkf_data_type;
@@ -50,21 +52,11 @@ typedef struct rml_enkf_data_struct rml_enkf_data_type;
set_subspace_dimension() routines will set one variable, AND
INVALIDATE THE OTHER. For most situations this will be OK, but if
you have repeated calls to both of these functions the end result
might be a surprise.
might be a surprise.
*/
#define INVALID_SUBSPACE_DIMENSION -1
#define INVALID_TRUNCATION -1
#define DEFAULT_SUBSPACE_DIMENSION INVALID_SUBSPACE_DIMENSION
#define DEFAULT_USE_PRIOR true
#define DEFAULT_LAMBDA_INCREASE_FACTOR 4
#define DEFAULT_LAMBDA_REDUCE_FACTOR 0.1
#define DEFAULT_LAMBDA0 -1
#define DEFAULT_LAMBDA_MIN 0.01
#define DEFAULT_LAMBDA_RECALCULATE false
#define DEFAULT_LOG_FILE "rml_enkf.out"
#define DEFAULT_CLEAR_LOG true
#define USE_PRIOR_KEY "USE_PRIOR"
#define LAMBDA_REDUCE_FACTOR_KEY "LAMBDA_REDUCE"
@@ -74,7 +66,7 @@ typedef struct rml_enkf_data_struct rml_enkf_data_type;
#define LAMBDA_RECALCULATE_KEY "LAMBDA_RECALCULATE"
#define ITER_KEY "ITER"
#define LOG_FILE_KEY "LOG_FILE"
#define CLEAR_LOG_KEY "CLEAR_LOG"
#define CLEAR_LOG_KEY "CLEAR_LOG"
@@ -103,32 +95,31 @@ typedef struct rml_enkf_data_struct rml_enkf_data_type;
destructor or free() function registered with the .freef field of
the analysis table.
*/
struct rml_enkf_data_struct {
UTIL_TYPE_ID_DECLARATION;
double truncation; // Controlled by config key: ENKF_TRUNCATION_KEY
int subspace_dimension; // Controlled by config key: ENKF_NCOMP_KEY (-1: use Truncation instead)
long option_flags;
int iteration_nr; // Keep track of the outer iteration loop
double Sk; // Objective function value
double Std; // Standard Deviation of the Objective function
double * Csc;
bool_vector_type * ens_mask;
double * Csc; // Vector with scalings for non-dimensionalizing states
matrix_type *Am; // Scaled right singular vectors of ensemble anomalies.
matrix_type *prior; // m_pr
matrix_type *state; // m_l
bool_vector_type * ens_mask; // Tells you which of the realisations are in use.
bool use_prior; // Use exact/approximate scheme? Approximate scheme drops the "prior" term in the LM step.
double lambda; // parameter to control the setp length in Marquardt levenberg optimization
double lambda0;
double lambda_min;
double lambda_reduce_factor;
double lambda_increase_factor;
bool lambda_recalculate;
bool clear_log;
char * log_file;
FILE * log_stream;
matrix_type *global_prior; // m_pr
matrix_type *previous_state; // m_l
double lambda; // parameter to control the setp length in Marquardt levenberg optimization
rml_enkf_log_type * rml_log;
rml_enkf_config_type * config;
};
@@ -144,69 +135,7 @@ static UTIL_SAFE_CAST_FUNCTION_CONST( rml_enkf_data , RML_ENKF_TYPE_ID )
//**********************************************
// Set / Get
//**********************************************
double rml_enkf_get_truncation( rml_enkf_data_type * data ) {
return data->truncation;
}
int rml_enkf_get_subspace_dimension( rml_enkf_data_type * data ) {
return data->subspace_dimension;
}
void rml_enkf_set_truncation( rml_enkf_data_type * data , double truncation ) {
data->truncation = truncation;
if (truncation > 0.0)
data->subspace_dimension = INVALID_SUBSPACE_DIMENSION;
}
void rml_enkf_set_lambda0( rml_enkf_data_type * data , double lambda0) {
data->lambda0 = lambda0;
}
double rml_enkf_get_lambda0( const rml_enkf_data_type * data ) {
return data->lambda0;
}
void rml_enkf_set_lambda_min( rml_enkf_data_type * data , double lambda_min) {
data->lambda_min = lambda_min;
}
double rml_enkf_get_lambda_min( const rml_enkf_data_type * data ) {
return data->lambda_min;
}
void rml_enkf_set_lambda_increase_factor( rml_enkf_data_type * data , double increase_factor) {
data->lambda_increase_factor = increase_factor;
}
double rml_enkf_get_lambda_increase_factor( const rml_enkf_data_type * data ) {
return data->lambda_increase_factor;
}
void rml_enkf_set_lambda_reduce_factor( rml_enkf_data_type * data , double reduce_factor) {
data->lambda_reduce_factor = reduce_factor;
}
double rml_enkf_get_lambda_reduce_factor( const rml_enkf_data_type * data ) {
return data->lambda_reduce_factor;
}
bool rml_enkf_get_use_prior( const rml_enkf_data_type * data ) {
return data->use_prior;
}
void rml_enkf_set_use_prior( rml_enkf_data_type * data , bool use_prior) {
data->use_prior = use_prior;
}
void rml_enkf_set_lambda_recalculate( rml_enkf_data_type * data , bool lambda_recalculate) {
data->lambda_recalculate = lambda_recalculate;
}
void rml_enkf_set_subspace_dimension( rml_enkf_data_type * data , int subspace_dimension) {
data->subspace_dimension = subspace_dimension;
if (subspace_dimension > 0)
data->truncation = INVALID_TRUNCATION;
}
void rml_enkf_set_iteration_nr( rml_enkf_data_type * data , int iteration_nr) {
data->iteration_nr = iteration_nr;
@@ -223,72 +152,36 @@ int rml_enkf_get_iteration_nr( const rml_enkf_data_type * data ) {
//**********************************************
// Log-file related stuff
//**********************************************
bool rml_enkf_get_clear_log( const rml_enkf_data_type * data ) {
return data->clear_log;
}
void rml_enkf_set_clear_log( rml_enkf_data_type * data , bool clear_log) {
data->clear_log = clear_log;
}
void rml_enkf_set_log_file( rml_enkf_data_type * data , const char * log_file ) {
data->log_file = util_realloc_string_copy( data->log_file , log_file );
}
const char * rml_enkf_get_log_file( const rml_enkf_data_type * data) {
return data->log_file;
}
void rml_enkf_log_line( rml_enkf_data_type * data , const char * fmt , ...) {
if (data->log_stream) {
va_list ap;
va_start(ap , fmt);
vfprintf( data->log_stream , fmt , ap );
va_end( ap );
}
}
static void rml_enkf_write_log_header( rml_enkf_data_type * data, const char * format) {
if (data->log_stream) {
if (rml_enkf_log_is_open( data->rml_log )) {
const char * column1 = "Iter#";
const char * column2 = "Lambda";
const char * column3 = "Sk old";
const char * column4 = "Sk_new";
const char * column5 = "std(Sk)";
rml_enkf_log_line(data, format, column1, column2, column3, column4, column5);
rml_enkf_log_line(data->rml_log, format, column1, column2, column3, column4, column5);
}
}
static void rml_enkf_write_iter_info( rml_enkf_data_type * data , double Sk_new, double Std_new ) {
if (data->log_stream) {
static void rml_enkf_write_iter_info( rml_enkf_data_type * data , double prev_Sk , double Sk_new, double Std_new ) {
if (rml_enkf_log_is_open( data->rml_log )) {
const char * format = "\n%2d-->%-2d %-7.3f %-7.3f --> %-7.3f %-7.3f";
const char * format_headers = "\n%-7s %-7s %-7s --> %-7s %-7s";
static int has_printed_header = 0;
static bool has_printed_header = false;
if (!has_printed_header) {
rml_enkf_write_log_header( data, format_headers );
has_printed_header = 1;
has_printed_header = true;
}
rml_enkf_log_line( data , format, data->iteration_nr, data->iteration_nr+1, data->lambda, data->Sk, Sk_new, Std_new);
rml_enkf_log_line( data->rml_log , format, data->iteration_nr, data->iteration_nr+1, data->lambda, prev_Sk, Sk_new, Std_new);
}
}
static void rml_enkf_open_log_file( rml_enkf_data_type * data ) {
data->log_stream = NULL;
if (data->log_file) {
if ( data->iteration_nr == 0) {
if (data->clear_log){
data->log_stream = util_mkdir_fopen( data->log_file , "w");
}
else
data->log_stream = util_mkdir_fopen( data->log_file , "a");
} else
data->log_stream = util_fopen( data->log_file , "a");
}
}
@@ -299,37 +192,28 @@ static void rml_enkf_open_log_file( rml_enkf_data_type * data ) {
void * rml_enkf_data_alloc( rng_type * rng) {
rml_enkf_data_type * data = util_malloc( sizeof * data);
UTIL_TYPE_ID_INIT( data , RML_ENKF_TYPE_ID );
data->log_file = NULL;
rml_enkf_set_truncation( data , DEFAULT_ENKF_TRUNCATION_ );
rml_enkf_set_subspace_dimension( data , DEFAULT_SUBSPACE_DIMENSION );
rml_enkf_set_use_prior( data , DEFAULT_USE_PRIOR );
rml_enkf_set_lambda0( data , DEFAULT_LAMBDA0 );
rml_enkf_set_lambda_increase_factor(data , DEFAULT_LAMBDA_INCREASE_FACTOR);
rml_enkf_set_lambda_reduce_factor(data , DEFAULT_LAMBDA_REDUCE_FACTOR);
rml_enkf_set_lambda_min( data , DEFAULT_LAMBDA_MIN );
rml_enkf_set_log_file( data , DEFAULT_LOG_FILE );
rml_enkf_set_clear_log( data , DEFAULT_CLEAR_LOG );
rml_enkf_set_lambda_recalculate( data , DEFAULT_LAMBDA_RECALCULATE );
data->config = rml_enkf_config_alloc();
data->rml_log = rml_enkf_log_alloc();
data->option_flags = ANALYSIS_NEED_ED + ANALYSIS_UPDATE_A + ANALYSIS_ITERABLE + ANALYSIS_SCALE_DATA;
data->Csc = NULL;
data->iteration_nr = 0;
data->Std = 0;
data->ens_mask = bool_vector_alloc(0,false);
data->state = matrix_alloc(1,1);
data->prior = matrix_alloc(1,1);
data->Std = 0;
data->previous_state = matrix_alloc(1,1);
data->global_prior = NULL;
data->ens_mask = NULL;
return data;
}
void rml_enkf_data_free( void * arg ) {
void rml_enkf_data_free( void * arg ) {
rml_enkf_data_type * data = rml_enkf_data_safe_cast( arg );
matrix_free( data->state );
matrix_free( data->prior );
matrix_free( data->previous_state );
if (data->global_prior)
matrix_free( data->global_prior );
util_safe_free( data->log_file );
bool_vector_free( data->ens_mask );
rml_enkf_log_free( data->rml_log );
rml_enkf_config_free( data->config );
free( data );
}
@@ -346,11 +230,10 @@ void rml_enkf_data_free( void * arg ) {
* Variable name in code <-> D.Oliver notation <-> Description
* -------------------------------------------------------------------------------------------------------------
* A <-> m_l <-> Ensemble matrix. Updated in-place by iterations.
* data->state <-> m_(l-1) <-> "A" from the previous iteration. Backs up A in case the update is bad.
* data->prior <-> <-> Previously: "active_prior". Stores A from before iter0, i.e. the actual prior.
* data->previous_state <-> m_(l-1) <-> "A" from the previous iteration. Backs up A in case the update is bad.
* data->global_prior <-> <-> Previously: "active_prior". Stores A from before iter0, i.e. the actual prior.
* Acopy <-> <-> Eliminated from code. Copy of A (at each iteration, before acceptance/rejection decision)
* data->prior0 <-> m_pr <-> Eliminated from code. Same as prior, but also includes columns (j) for which ens_mask[j]==false.
* Seems pointless. Only creates confusion.
*
* Am <-> A_m <-> Am = Um*Wm^(-1)
* Csc <-> C_sc^(1/2) <-> State scalings. Note the square root.
@@ -373,28 +256,30 @@ void rml_enkf_data_free( void * arg ) {
static void rml_enkf_init1__( rml_enkf_data_type * data) {
// Differentiate this routine from init2__, which actually calculates the prior mismatch update.
// This routine does not change any ensemble matrix.
// Um*Wm^(-1) are the scaled, truncated, right singular vectors of data->prior
// Um*Wm^(-1) are the scaled, truncated, right singular vectors of data->global_prior
int state_size = matrix_get_rows( data->prior );
int ens_size = matrix_get_columns( data->prior );
int nrmin = util_int_min( ens_size , state_size);
matrix_type * Dm = matrix_alloc_copy( data->prior );
matrix_type * Um = matrix_alloc( state_size , nrmin ); /* Left singular vectors. */
matrix_type * VmT = matrix_alloc( nrmin , ens_size ); /* Right singular vectors. */
double * Wm = util_calloc( nrmin , sizeof * Wm );
double nsc = 1/sqrt(ens_size - 1);
matrix_type * prior = matrix_alloc_column_compressed_copy( data->global_prior , data->ens_mask);
int state_size = matrix_get_rows( prior );
int ens_size = matrix_get_columns( prior );
int nrmin = util_int_min( ens_size , state_size);
matrix_type * Dm = matrix_alloc_copy( prior );
matrix_type * Um = matrix_alloc( state_size , nrmin ); /* Left singular vectors. */
matrix_type * VmT = matrix_alloc( nrmin , ens_size ); /* Right singular vectors. */
double * Wm = util_calloc( nrmin , sizeof * Wm );
double nsc = 1/sqrt(ens_size - 1);
matrix_subtract_row_mean(Dm);
for (int i=0; i < state_size; i++){
double sc = nsc / (data->Csc[i]);
matrix_scale_row( Dm , i , sc);
{
const double * Csc = data->Csc;
for (int i=0; i < state_size; i++){
double sc = nsc / (Csc[i]);
matrix_scale_row( Dm , i , sc);
}
}
// Um Wm VmT = Dm; nsign1 = num of non-zero singular values.
int nsign1 = enkf_linalg_svd_truncation(Dm , data->truncation , -1 , DGESVD_MIN_RETURN , Wm , Um , VmT);
int nsign1 = enkf_linalg_svd_truncation(Dm , rml_enkf_config_get_truncation( data->config ) , -1 , DGESVD_MIN_RETURN , Wm , Um , VmT);
// Am = Um*Wm^(-1). I.e. scale *columns* of Um
enkf_linalg_rml_enkfAm(Um, Wm, nsign1);
@@ -402,25 +287,31 @@ static void rml_enkf_init1__( rml_enkf_data_type * data) {
matrix_free(Um);
matrix_free(VmT);
matrix_free(Dm);
matrix_free(prior);
free(Wm);
}
// Creates state scaling matrix
void rml_enkf_init_Csc(rml_enkf_data_type * data){
void rml_enkf_init_Csc(const rml_enkf_data_type * data ){
// This seems a strange choice of scaling matrix. Review?
int state_size = matrix_get_rows( data->prior );
int ens_size = matrix_get_columns( data->prior );
matrix_type * prior = matrix_alloc_column_compressed_copy( data->global_prior , data->ens_mask );
{
int state_size = matrix_get_rows( prior );
int ens_size = matrix_get_columns( prior );
for (int row=0; row < state_size; row++) {
double sumrow = matrix_get_row_sum(data->prior , row);
double tmp = sumrow / ens_size;
for (int row=0; row < state_size; row++) {
double sumrow = matrix_get_row_sum(prior , row);
double tmp = sumrow / ens_size;
if (abs(tmp)< 1)
data->Csc[row] = 0.05;
else
data->Csc[row] = 1.00;
if (abs(tmp)< 1)
data->Csc[row] = 0.05;
else
data->Csc[row] = 1.00;
}
matrix_free( prior );
}
}
@@ -428,6 +319,7 @@ void rml_enkf_init_Csc(rml_enkf_data_type * data){
static void rml_enkf_initA__(rml_enkf_data_type * data, matrix_type * A, matrix_type * S, matrix_type * Cd, matrix_type * E, matrix_type * D, matrix_type * Udr, double * Wdr, matrix_type * VdTr) {
int ens_size = matrix_get_columns( S );
int state_size = matrix_get_rows( A );
double nsc = 1/sqrt(ens_size-1);
int nsign;
@@ -440,37 +332,37 @@ static void rml_enkf_initA__(rml_enkf_data_type * data, matrix_type * A, matrix_
matrix_matmul(tmp , Cd , S ); //
matrix_scale(tmp , nsc); //
nsign = enkf_linalg_svd_truncation(tmp , data->truncation , -1 , DGESVD_MIN_RETURN , Wdr , Udr , VdTr);
nsign = enkf_linalg_svd_truncation(tmp , rml_enkf_config_get_truncation( data->config ) , -1 , DGESVD_MIN_RETURN , Wdr , Udr , VdTr);
matrix_free( tmp );
}
// Calc X3
{
matrix_type * X3 = matrix_alloc( ens_size, ens_size );
{
matrix_type * X1 = matrix_alloc( nsign, ens_size);
matrix_type * X2 = matrix_alloc( nsign, ens_size );
// See LM-EnRML algorithm in Oliver'2013 (Comp. Geo.) for meaning
enkf_linalg_rml_enkfX1(X1, Udr ,D ,Cd ); // X1 = Ud(T)*Cd(-1/2)*D -- D= -(dk-d0)
enkf_linalg_rml_enkfX2(X2, Wdr ,X1 ,data->lambda + 1 , nsign); // X2 = ((a*Ipd)+Wd^2)^-1 * X1
enkf_linalg_rml_enkfX3(X3, VdTr ,Wdr,X2, nsign); // X3 = Vd *Wd*X2
matrix_free(X2);
matrix_free(X1);
}
// Update A
{
matrix_type * dA1 = matrix_alloc( matrix_get_rows(A) , ens_size);
matrix_type * Dm = matrix_alloc_copy( A );
matrix_type * dA1 = matrix_alloc( state_size , ens_size);
matrix_type * Dm = matrix_alloc_copy( A );
matrix_subtract_row_mean( Dm ); /* Remove the mean from the ensemble of model parameters*/
matrix_scale(Dm, nsc);
matrix_matmul(dA1, Dm , X3);
matrix_inplace_add(A,dA1); // dA
matrix_inplace_add(A,dA1); // dA
matrix_free(Dm);
matrix_free(dA1);
@@ -487,15 +379,15 @@ void rml_enkf_init2__( rml_enkf_data_type * data, matrix_type *A, double * Wdr,
int state_size = matrix_get_rows( A );
int ens_size = matrix_get_columns( A );
double nsc = 1/sqrt(ens_size-1);
double nsc = 1/sqrt(ens_size-1);
matrix_type *Am = matrix_alloc_copy(data->Am);
matrix_type *Apr = matrix_alloc_copy(data->prior);
matrix_type *Apr = matrix_alloc_column_compressed_copy(data->global_prior , data->ens_mask );
// fprintf(stdout,"\n");
// fprintf(stdout,"A: %d x %d\n", matrix_get_rows(A), matrix_get_columns(A));
// fprintf(stdout,"prior : %d x %d\n", matrix_get_rows(data->prior), matrix_get_columns(data->prior));
// fprintf(stdout,"state : %d x %d\n", matrix_get_rows(data->state), matrix_get_columns(data->state));
// fprintf(stdout,"prior : %d x %d\n", matrix_get_rows(data->global_prior), matrix_get_columns(data->global_prior));
// fprintf(stdout,"state : %d x %d\n", matrix_get_rows(data->previous_state), matrix_get_columns(data->previous_state));
// fprintf(stdout,"Apr : %d x %d\n", matrix_get_rows(Apr), matrix_get_columns(Apr));
// fprintf(stdout,"Am : %d x %d\n", matrix_get_rows(Am), matrix_get_columns(Am));
// Example:
@@ -508,7 +400,7 @@ void rml_enkf_init2__( rml_enkf_data_type * data, matrix_type *A, double * Wdr,
int nsign1 = matrix_get_columns(data->Am);
matrix_type * X4 = matrix_alloc(nsign1,ens_size);
matrix_type * X5 = matrix_alloc(state_size,ens_size);
@@ -516,13 +408,15 @@ void rml_enkf_init2__( rml_enkf_data_type * data, matrix_type *A, double * Wdr,
matrix_type * X7 = matrix_alloc(ens_size,ens_size);
matrix_type * dA2 = matrix_alloc(state_size , ens_size);
matrix_type * Dk1 = matrix_alloc_copy( A );
// Dk = Csc^(-1) * (A - Aprior)
// X4 = Am' * Dk
{
matrix_type * Dk = matrix_alloc_copy( A );
matrix_inplace_sub(Dk, Apr);
matrix_inplace_sub( Dk , Apr );
rml_enkf_common_scaleA(Dk , data->Csc , true);
matrix_dgemm(X4 , Am , Dk , true, false, 1.0, 0.0);
matrix_free(Dk);
}
@@ -536,10 +430,10 @@ void rml_enkf_init2__( rml_enkf_data_type * data, matrix_type *A, double * Wdr,
// X6 = Dk1' * X5
matrix_dgemm(X6, Dk1, X5, true, false, 1.0, 0.0);
// X7
enkf_linalg_rml_enkfX7(X7, VdTr , Wdr , data->lambda + 1, X6);
// delta m_2
rml_enkf_common_scaleA(Dk1 , data->Csc , false);
matrix_matmul(dA2 , Dk1 , X7);
@@ -547,7 +441,7 @@ void rml_enkf_init2__( rml_enkf_data_type * data, matrix_type *A, double * Wdr,
matrix_free(Am);
matrix_free(Apr);
matrix_free(X4);
matrix_free(X4);
matrix_free(X5);
matrix_free(X6);
matrix_free(X7);
@@ -557,39 +451,45 @@ void rml_enkf_init2__( rml_enkf_data_type * data, matrix_type *A, double * Wdr,
// Initialize state and prior from A. Initialize lambda0, lambda. Call initA__, init1__
static void rml_enkf_updateA_iter0(rml_enkf_data_type * data, matrix_type * A, matrix_type * S, matrix_type * R, matrix_type * dObs, matrix_type * E, matrix_type * D, matrix_type * Cd) {
int ens_size = matrix_get_columns( S );
int nrobs = matrix_get_rows( S );
int nrmin = util_int_min( ens_size , nrobs);
int nrmin = util_int_min( ens_size , nrobs);
int state_size = matrix_get_rows( A );
matrix_type * Skm = matrix_alloc(ens_size, ens_size); // Mismatch
matrix_type * Ud = matrix_alloc( nrobs , nrmin ); /* Left singular vectors. */
matrix_type * VdT = matrix_alloc( nrmin , ens_size ); /* Right singular vectors. */
double * Wd = util_calloc( nrmin , sizeof * Wd );
double * Wd = util_calloc( nrmin , sizeof * Wd );
data->Csc = util_calloc(state_size , sizeof * data->Csc);
data->Sk = enkf_linalg_data_mismatch(D,Cd,Skm);
data->Sk = enkf_linalg_data_mismatch(D,Cd,Skm);
data->Std = matrix_diag_std(Skm,data->Sk);
if (data->lambda0 < 0)
data->lambda = pow(10 , floor(log10(data->Sk/(2*nrobs))) );
else
data->lambda = data->lambda0;
// state = A, prior = A
rml_enkf_common_store_state( data->state , A , data->ens_mask );
rml_enkf_common_recover_state( A , data->prior , data->ens_mask );
{
double lambda0 = rml_enkf_config_get_lambda0( data->config );
if (lambda0 < 0)
data->lambda = pow(10 , floor(log10(data->Sk/(2*nrobs))) );
else
data->lambda = lambda0;
}
// state = A
rml_enkf_common_store_state( data->previous_state , A , data->ens_mask );
// prior = A
data->global_prior = matrix_alloc_copy( data->previous_state );
// Update dependant on data mismatch
rml_enkf_initA__(data , A, S , Cd , E , D , Ud , Wd , VdT);
// Update dependant on prior mismatch. This should be zero (coz iter0).
// Therefore the purpose of init1__ is just to prepare some matrices.
if (data->use_prior) {
if (rml_enkf_config_get_use_prior(data->config)) {
rml_enkf_init_Csc( data );
rml_enkf_init1__(data );
rml_enkf_init1__( data );
}
rml_enkf_write_iter_info(data, data->Sk, data->Std);
rml_enkf_write_iter_info(data, data->Sk , data->Sk, data->Std);
matrix_free( Skm );
matrix_free( Ud );
@@ -597,7 +497,7 @@ static void rml_enkf_updateA_iter0(rml_enkf_data_type * data, matrix_type * A, m
free( Wd );
}
// Main routine. Controls the iterations. Called from analysis_module.c: analysis_module_updateA()
void rml_enkf_updateA(void * module_data, matrix_type * A, matrix_type * S, matrix_type * R, matrix_type * dObs, matrix_type * E, matrix_type * D) {
// A : ensemble matrix
// R : (Inv?) Obs error cov.
@@ -616,12 +516,12 @@ void rml_enkf_updateA(void * module_data, matrix_type * A, matrix_type * S, matr
double nsc = 1/sqrt(ens_size-1); // Scale factor
matrix_type * Cd = matrix_alloc( nrobs, nrobs ); // Cov(E), where E = measurement perturbations?
// Empirical error covar. R is left unused. Investigate?
// Empirical error covar. R is left unused. Investigate?
enkf_linalg_Covariance(Cd ,E ,nsc, nrobs); // Cd = SampCov(E) (including (N-1) normalization)
matrix_inv(Cd); // In-place inversion
rml_enkf_open_log_file(data);
rml_enkf_log_open(data->rml_log , data->iteration_nr);
fprintf(stdout,"\nIter %d --> %d", data->iteration_nr, data->iteration_nr + 1);
@@ -641,9 +541,9 @@ void rml_enkf_updateA(void * module_data, matrix_type * A, matrix_type * S, matr
// Lambda = Normalized data mismatch (rounded)
if (data->lambda_recalculate)
if (rml_enkf_config_get_lambda_recalculate( data->config ))
data->lambda = pow(10 , floor(log10(Sk_new / (2*nrobs))) );
// Accept/Reject update? Lambda calculation.
{
bool mismatch_reduced = false;
@@ -651,65 +551,70 @@ void rml_enkf_updateA(void * module_data, matrix_type * A, matrix_type * S, matr
if (Sk_new < data->Sk)
mismatch_reduced = true;
if (Std_new <= data->Std)
std_reduced = true;
fprintf(stdout,"\nWriting iter info to file now. Iter %d --> %d", data->iteration_nr, data->iteration_nr + 1);
rml_enkf_write_iter_info(data, Sk_new, Std_new);
rml_enkf_write_iter_info(data, data->Sk , Sk_new, Std_new);
if (mismatch_reduced) {
/*
Stop check: if ( (1- (Sk_new/data->Sk)) < .0001) // check convergence ** model change norm has to be added in this!!
*/
// Reduce Lambda
if (std_reduced)
data->lambda = data->lambda * data->lambda_reduce_factor;
rml_enkf_common_store_state(data->state , A , data->ens_mask );
if (std_reduced)
data->lambda = data->lambda * rml_enkf_config_get_lambda_decrease_factor( data->config );
rml_enkf_common_store_state(data->previous_state , A , data->ens_mask );
data->Sk = Sk_new;
data->Std=Std_new;
data->iteration_nr++;
} else {
// Increase lambda
data->lambda = data->lambda * data->lambda_increase_factor;
// A = data->state
rml_enkf_common_recover_state( data->state , A , data->ens_mask );
data->lambda = data->lambda * rml_enkf_config_get_lambda_increase_factor( data->config );
// A = data->previous_state
rml_enkf_common_recover_state( data->previous_state , A , data->ens_mask );
}
}
// Update dependant on data mismatch (delta m_1)
rml_enkf_initA__(data , A , S , Cd , E , D , Ud , Wd , VdT);
// Update dependant on prior mismatch (delta m_2)
if (data->use_prior) {
if (rml_enkf_config_get_use_prior(data->config)) {
rml_enkf_init_Csc( data );
rml_enkf_init2__(data , A , Wd , VdT);
}
// Free
matrix_free(Skm);
matrix_free( Ud );
matrix_free( VdT );
free( Wd );
}
if (data->lambda < data->lambda_min)
data->lambda = data->lambda_min;
{
double lambda_min = rml_enkf_config_get_lambda_min( data->config );
if (data->lambda < lambda_min)
data->lambda = lambda_min;
}
if (data->log_stream)
fclose( data->log_stream );
rml_enkf_log_close( data->rml_log );
matrix_free(Cd);
}
// Called from analysis_module.c: analysis_module_init_update()
void rml_enkf_init_update(void * arg, const bool_vector_type * ens_mask, const matrix_type * S, const matrix_type * R, const matrix_type * dObs, const matrix_type * E, const matrix_type * D ) {
void rml_enkf_init_update(void * arg, const bool_vector_type * ens_mask, const matrix_type * S, const matrix_type * R, const matrix_type * dObs, const matrix_type * E, const matrix_type * D ) {
rml_enkf_data_type * module_data = rml_enkf_data_safe_cast( arg );
bool_vector_memcpy( module_data->ens_mask , ens_mask );
if (module_data->ens_mask)
bool_vector_free( module_data->ens_mask );
module_data->ens_mask = bool_vector_alloc_copy( ens_mask );
}
@@ -717,6 +622,7 @@ void rml_enkf_init_update(void * arg, const bool_vector_type * ens_mask, const m
//**********************************************
// Set / Get basic types
//**********************************************
@@ -724,9 +630,9 @@ bool rml_enkf_set_int( void * arg , const char * var_name , int value) {
rml_enkf_data_type * module_data = rml_enkf_data_safe_cast( arg );
{
bool name_recognized = true;
if (strcmp( var_name , ENKF_NCOMP_KEY_) == 0)
rml_enkf_set_subspace_dimension( module_data , value );
rml_enkf_config_set_subspace_dimension(module_data->config , value);
else if (strcmp( var_name , ITER_KEY) == 0)
rml_enkf_set_iteration_nr( module_data , value );
else
@@ -750,13 +656,13 @@ bool rml_enkf_set_bool( void * arg , const char * var_name , bool value) {
rml_enkf_data_type * module_data = rml_enkf_data_safe_cast( arg );
{
bool name_recognized = true;
if (strcmp( var_name , USE_PRIOR_KEY) == 0)
rml_enkf_set_use_prior( module_data , value );
rml_enkf_config_set_use_prior( module_data->config , value);
else if (strcmp( var_name , CLEAR_LOG_KEY) == 0)
rml_enkf_set_clear_log( module_data , value );
rml_enkf_log_set_clear_log( module_data->rml_log , value );
else if (strcmp( var_name , LAMBDA_RECALCULATE_KEY) == 0)
rml_enkf_set_lambda_recalculate( module_data , value );
rml_enkf_config_set_lambda_recalculate( module_data->config , value );
else
name_recognized = false;
@@ -768,11 +674,11 @@ bool rml_enkf_get_bool( const void * arg, const char * var_name) {
const rml_enkf_data_type * module_data = rml_enkf_data_safe_cast_const( arg );
{
if (strcmp(var_name , USE_PRIOR_KEY) == 0)
return module_data->use_prior;
else if (strcmp(var_name , CLEAR_LOG_KEY) == 0)
return module_data->clear_log;
else if (strcmp(var_name , LAMBDA_RECALCULATE_KEY) == 0)
return module_data->lambda_recalculate;
return rml_enkf_config_get_use_prior( module_data->config );
else if (strcmp(var_name , CLEAR_LOG_KEY) == 0)
return rml_enkf_log_get_clear_log( module_data->rml_log );
else if (strcmp(var_name , LAMBDA_RECALCULATE_KEY) == 0)
return rml_enkf_config_get_lambda_recalculate( module_data->config );
else
return false;
}
@@ -784,15 +690,15 @@ bool rml_enkf_set_double( void * arg , const char * var_name , double value) {
bool name_recognized = true;
if (strcmp( var_name , ENKF_TRUNCATION_KEY_) == 0)
rml_enkf_set_truncation( module_data , value );
rml_enkf_config_set_truncation( module_data->config , value );
else if (strcmp( var_name , LAMBDA_INCREASE_FACTOR_KEY) == 0)
rml_enkf_set_lambda_increase_factor( module_data , value );
rml_enkf_config_set_lambda_increase_factor( module_data->config , value );
else if (strcmp( var_name , LAMBDA_REDUCE_FACTOR_KEY) == 0)
rml_enkf_set_lambda_reduce_factor( module_data , value );
rml_enkf_config_set_lambda_decrease_factor( module_data->config , value );
else if (strcmp( var_name , LAMBDA0_KEY) == 0)
rml_enkf_set_lambda0( module_data , value );
rml_enkf_config_set_lambda0( module_data->config , value );
else if (strcmp( var_name , LAMBDA_MIN_KEY) == 0)
rml_enkf_set_lambda_min( module_data , value );
rml_enkf_config_set_lambda_min( module_data->config , value );
else
name_recognized = false;
@@ -804,17 +710,21 @@ double rml_enkf_get_double( const void * arg, const char * var_name) {
const rml_enkf_data_type * module_data = rml_enkf_data_safe_cast_const( arg );
{
if (strcmp(var_name , LAMBDA_REDUCE_FACTOR_KEY) == 0)
return module_data->lambda_reduce_factor;
return rml_enkf_config_get_lambda_decrease_factor(module_data->config);
if (strcmp(var_name , LAMBDA_INCREASE_FACTOR_KEY) == 0)
return module_data->lambda_increase_factor;
return rml_enkf_config_get_lambda_increase_factor(module_data->config);
if (strcmp(var_name , LAMBDA0_KEY) == 0)
return module_data->lambda0;
return rml_enkf_config_get_lambda0(module_data->config);
if (strcmp(var_name , LAMBDA_MIN_KEY) == 0)
return module_data->lambda_min;
return rml_enkf_config_get_lambda_min(module_data->config);
if (strcmp(var_name , ENKF_TRUNCATION_KEY_) == 0)
return module_data->truncation;
else
return -1;
return rml_enkf_config_get_truncation( module_data->config );
return -1;
}
}
@@ -823,9 +733,9 @@ bool rml_enkf_set_string( void * arg , const char * var_name , const char * valu
rml_enkf_data_type * module_data = rml_enkf_data_safe_cast( arg );
{
bool name_recognized = true;
if (strcmp( var_name , LOG_FILE_KEY) == 0)
rml_enkf_set_log_file( module_data , value );
rml_enkf_log_set_log_file( module_data->rml_log , value );
else
name_recognized = false;
@@ -836,7 +746,7 @@ bool rml_enkf_set_string( void * arg , const char * var_name , const char * valu
long rml_enkf_get_options( void * arg , long flag ) {
rml_enkf_data_type * module_data = rml_enkf_data_safe_cast( arg );
{
return module_data->option_flags;
return rml_enkf_config_get_option_flags( module_data->config );
}
}
@@ -871,7 +781,7 @@ void * rml_enkf_get_ptr( const void * arg , const char * var_name ) {
const rml_enkf_data_type * module_data = rml_enkf_data_safe_cast_const( arg );
{
if (strcmp(var_name , LOG_FILE_KEY) == 0)
return module_data->log_file;
return (void *) rml_enkf_log_get_log_file( module_data->rml_log );
else
return NULL;
}
@@ -895,13 +805,13 @@ void * rml_enkf_get_ptr( const void * arg , const char * var_name ) {
analysis_table_type SYMBOL_TABLE = {
.alloc = rml_enkf_data_alloc,
.freef = rml_enkf_data_free,
.set_int = rml_enkf_set_int ,
.set_double = rml_enkf_set_double ,
.set_bool = rml_enkf_set_bool,
.set_int = rml_enkf_set_int ,
.set_double = rml_enkf_set_double ,
.set_bool = rml_enkf_set_bool,
.set_string = rml_enkf_set_string,
.get_options = rml_enkf_get_options ,
.get_options = rml_enkf_get_options ,
.initX = NULL,
.updateA = rml_enkf_updateA ,
.updateA = rml_enkf_updateA ,
.init_update = rml_enkf_init_update ,
.complete_update = NULL,
.has_var = rml_enkf_has_var,

36
ThirdParty/Ert/devel/libanalysis/modules/rml_enkf_common.c vendored
View File

@@ -1,19 +1,19 @@
/*
Copyright (C) 2011 Statoil ASA, Norway.
The file 'rml_enkf.c' 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.
Copyright (C) 2011 Statoil ASA, Norway.
The file 'rml_enkf.c' 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.
*/
#include <stdlib.h>
@@ -38,7 +38,7 @@
// zzz_enkf_common_recover_state(state , A ,ens_mask) assigns state to A. RESIZES A to rows(state)-by-SUM(ens_mask)
void rml_enkf_common_store_state( matrix_type * state , const matrix_type * A , const bool_vector_type * ens_mask ) {
void rml_enkf_common_store_state( matrix_type * state , const matrix_type * A , const bool_vector_type * ens_mask ) {
matrix_resize( state , matrix_get_rows( A ) , bool_vector_size( ens_mask ) , false);
{
const int ens_size = bool_vector_size( ens_mask );
@@ -55,11 +55,11 @@ void rml_enkf_common_store_state( matrix_type * state , const matrix_type * A ,
void rml_enkf_common_recover_state( const matrix_type * state , matrix_type * A , const bool_vector_type * ens_mask ) {
void rml_enkf_common_recover_state( const matrix_type * state , matrix_type * A , const bool_vector_type * ens_mask ) {
const int ens_size = bool_vector_size( ens_mask );
const int active_size = bool_vector_count_equal( ens_mask , true );
const int rows = matrix_get_rows( state );
matrix_resize( A , rows , active_size , false );
{
int active_index = 0;

168
ThirdParty/Ert/devel/libanalysis/modules/rml_enkf_config.c vendored Normal file
View File

@@ -0,0 +1,168 @@
/*
Copyright (C) 2015 Statoil ASA, Norway.
The file 'rml_enkf_config.c' 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.
*/
#include <ert/util/util.h>
#include <ert/util/type_macros.h>
#include <ert/analysis/std_enkf.h>
#include <ert/analysis/analysis_module.h>
#include <rml_enkf_config.h>
#define INVALID_SUBSPACE_DIMENSION -1
#define INVALID_TRUNCATION -1
#define DEFAULT_SUBSPACE_DIMENSION INVALID_SUBSPACE_DIMENSION
#define DEFAULT_USE_PRIOR true
#define DEFAULT_LAMBDA_INCREASE_FACTOR 4
#define DEFAULT_LAMBDA_REDUCE_FACTOR 0.1
#define DEFAULT_LAMBDA0 -1
#define DEFAULT_LAMBDA_MIN 0.01
#define DEFAULT_LAMBDA_RECALCULATE false
#define RML_ENKF_CONFIG_TYPE_ID 61400061
struct rml_enkf_config_struct {
UTIL_TYPE_ID_DECLARATION;
double truncation; // Controlled by config key: ENKF_TRUNCATION_KEY
int subspace_dimension; // Controlled by config key: ENKF_NCOMP_KEY (-1: use Truncation instead)
long option_flags;
bool use_prior; // Use exact/approximate scheme? Approximate scheme drops the "prior" term in the LM step.
double lambda0;
double lambda_min;
double lambda_decrease_factor;
double lambda_increase_factor;
bool lambda_recalculate;
};
rml_enkf_config_type * rml_enkf_config_alloc() {
rml_enkf_config_type * config = util_malloc( sizeof * config );
UTIL_TYPE_ID_INIT( config , RML_ENKF_CONFIG_TYPE_ID );
rml_enkf_config_set_truncation( config , DEFAULT_ENKF_TRUNCATION_);
rml_enkf_config_set_subspace_dimension( config , DEFAULT_SUBSPACE_DIMENSION);
rml_enkf_config_set_use_prior( config , DEFAULT_USE_PRIOR );
rml_enkf_config_set_option_flags( config , ANALYSIS_NEED_ED + ANALYSIS_UPDATE_A + ANALYSIS_ITERABLE + ANALYSIS_SCALE_DATA);
rml_enkf_config_set_lambda_min( config , DEFAULT_LAMBDA_MIN );
rml_enkf_config_set_lambda0( config , DEFAULT_LAMBDA0 );
rml_enkf_config_set_lambda_decrease_factor( config , DEFAULT_LAMBDA_REDUCE_FACTOR );
rml_enkf_config_set_lambda_increase_factor( config , DEFAULT_LAMBDA_INCREASE_FACTOR );
rml_enkf_config_set_lambda_recalculate( config , DEFAULT_LAMBDA_RECALCULATE );
return config;
}
bool rml_enkf_config_get_use_prior( const rml_enkf_config_type * config ) {
return config->use_prior;
}
void rml_enkf_config_set_use_prior( rml_enkf_config_type * config , bool use_prior) {
config->use_prior = use_prior;
}
double rml_enkf_config_get_truncation( rml_enkf_config_type * config ) {
return config->truncation;
}
void rml_enkf_config_set_truncation( rml_enkf_config_type * config , double truncation) {
config->truncation = truncation;
if (truncation > 0.0)
config->subspace_dimension = INVALID_SUBSPACE_DIMENSION;
}
int rml_enkf_config_get_subspace_dimension( rml_enkf_config_type * config ) {
return config->subspace_dimension;
}
void rml_enkf_config_set_subspace_dimension( rml_enkf_config_type * config , int subspace_dimension) {
config->subspace_dimension = subspace_dimension;
if (subspace_dimension > 0)
config->truncation = INVALID_TRUNCATION;
}
void rml_enkf_config_set_option_flags( rml_enkf_config_type * config , long flags) {
config->option_flags = flags;
}
long rml_enkf_config_get_option_flags( const rml_enkf_config_type * config ) {
return config->option_flags;
}
double rml_enkf_config_get_lambda0( rml_enkf_config_type * config ) {
return config->lambda0;
}
void rml_enkf_config_set_lambda0( rml_enkf_config_type * config , double lambda0) {
config->lambda0 = lambda0;
}
double rml_enkf_config_get_lambda_min( rml_enkf_config_type * config ) {
return config->lambda_min;
}
void rml_enkf_config_set_lambda_min( rml_enkf_config_type * config , double lambda_min) {
config->lambda_min = lambda_min;
}
double rml_enkf_config_get_lambda_increase_factor( rml_enkf_config_type * config ) {
return config->lambda_increase_factor;
}
void rml_enkf_config_set_lambda_increase_factor( rml_enkf_config_type * config , double lambda_increase_factor) {
config->lambda_increase_factor = lambda_increase_factor;
}
double rml_enkf_config_get_lambda_decrease_factor( rml_enkf_config_type * config ) {
return config->lambda_decrease_factor;
}
void rml_enkf_config_set_lambda_decrease_factor( rml_enkf_config_type * config , double lambda_decrease_factor) {
config->lambda_decrease_factor = lambda_decrease_factor;
}
bool rml_enkf_config_get_lambda_recalculate( const rml_enkf_config_type * config ) {
return config->lambda_recalculate;
}
void rml_enkf_config_set_lambda_recalculate( rml_enkf_config_type * config , bool lambda_recalculate) {
config->lambda_recalculate = lambda_recalculate;
}
void rml_enkf_config_free(rml_enkf_config_type * config) {
free( config );
}

64
ThirdParty/Ert/devel/libanalysis/modules/rml_enkf_config.h vendored Normal file
View File

@@ -0,0 +1,64 @@
/*
Copyright (C) 2015 Statoil ASA, Norway.
The file 'rml_enkf_config.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 RML_ENKF_CONFIG_H
#define RML_ENKF_CONFIG_H
#ifdef __cplusplus
extern "C" {
#endif
typedef struct rml_enkf_config_struct rml_enkf_config_type;
rml_enkf_config_type * rml_enkf_config_alloc();
void rml_enkf_config_free(rml_enkf_config_type * config);
int rml_enkf_config_get_subspace_dimension( rml_enkf_config_type * config );
void rml_enkf_config_set_subspace_dimension( rml_enkf_config_type * config , int subspace_dimension);
double rml_enkf_config_get_truncation( rml_enkf_config_type * config );
void rml_enkf_config_set_truncation( rml_enkf_config_type * config , double truncation);
bool rml_enkf_config_get_use_prior( const rml_enkf_config_type * config );
void rml_enkf_config_set_use_prior( rml_enkf_config_type * config , bool use_prior);
void rml_enkf_config_set_option_flags( rml_enkf_config_type * config , long flags);
long rml_enkf_config_get_option_flags( const rml_enkf_config_type * config );
double rml_enkf_config_get_lambda0( rml_enkf_config_type * config );
void rml_enkf_config_set_lambda0( rml_enkf_config_type * config , double lambda0);
double rml_enkf_config_get_lambda_min( rml_enkf_config_type * config );
void rml_enkf_config_set_lambda_min( rml_enkf_config_type * config , double lambda_min);
double rml_enkf_config_get_lambda_increase_factor( rml_enkf_config_type * config );
void rml_enkf_config_set_lambda_increase_factor( rml_enkf_config_type * config , double lambda_increase_factor);
double rml_enkf_config_get_lambda_decrease_factor( rml_enkf_config_type * config );
void rml_enkf_config_set_lambda_decrease_factor( rml_enkf_config_type * config , double lambda_decrease_factor);
bool rml_enkf_config_get_lambda_recalculate( const rml_enkf_config_type * config );
void rml_enkf_config_set_lambda_recalculate( rml_enkf_config_type * config , bool lambda_recalculate);
#ifdef __cplusplus
}
#endif
#endif

108
ThirdParty/Ert/devel/libanalysis/modules/rml_enkf_log.c vendored Normal file
View File

@@ -0,0 +1,108 @@
/*
Copyright (C) 2015 Statoil ASA, Norway.
The file 'rml_enkf_log.c' 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.
*/
#include <stdio.h>
#include <stdlib.h>
#include <ert/util/util.h>
#include <rml_enkf_log.h>
#define DEFAULT_LOG_FILE "rml_enkf.out"
#define DEFAULT_CLEAR_LOG true
struct rml_enkf_log_struct {
bool clear_log;
char * log_file;
FILE * log_stream;
};
rml_enkf_log_type * rml_enkf_log_alloc() {
rml_enkf_log_type * rml_log = util_malloc( sizeof * rml_log );
rml_log->log_file = NULL;
rml_log->log_stream = NULL;
rml_enkf_log_set_clear_log( rml_log , DEFAULT_CLEAR_LOG );
return rml_log;
}
bool rml_enkf_log_get_clear_log( const rml_enkf_log_type * data ) {
return data->clear_log;
}
void rml_enkf_log_set_clear_log( rml_enkf_log_type * data , bool clear_log) {
data->clear_log = clear_log;
}
void rml_enkf_log_set_log_file( rml_enkf_log_type * data , const char * log_file ) {
data->log_file = util_realloc_string_copy( data->log_file , log_file );
}
const char * rml_enkf_log_get_log_file( const rml_enkf_log_type * data) {
return data->log_file;
}
void rml_enkf_log_free(rml_enkf_log_type * rml_log) {
rml_enkf_log_close( rml_log );
util_safe_free( rml_log->log_file );
free( rml_log );
}
void rml_enkf_log_open( rml_enkf_log_type * rml_log , int iteration_nr ) {
if (rml_log->log_file) {
if ( iteration_nr == 0) {
if (rml_log->clear_log)
rml_log->log_stream = util_mkdir_fopen( rml_log->log_file , "w");
else
rml_log->log_stream = util_mkdir_fopen( rml_log->log_file , "a");
} else
rml_log->log_stream = util_fopen( rml_log->log_file , "a");
}
}
bool rml_enkf_log_is_open( const rml_enkf_log_type * rml_log ) {
if (rml_log->log_stream)
return true;
else
return false;
}
void rml_enkf_log_close( rml_enkf_log_type * rml_log ) {
if (rml_log->log_stream)
fclose( rml_log->log_stream );
rml_log->log_stream = NULL;
}
void rml_enkf_log_line( rml_enkf_log_type * rml_log , const char * fmt , ...) {
if (rml_log->log_stream) {
va_list ap;
va_start(ap , fmt);
vfprintf( rml_log->log_stream , fmt , ap );
va_end( ap );
}
}

45
ThirdParty/Ert/devel/libanalysis/modules/rml_enkf_log.h vendored Normal file
View File

@@ -0,0 +1,45 @@
/*
Copyright (C) 2015 Statoil ASA, Norway.
The file 'rml_enkf_log.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 RML_ENKF_LOG_H
#define RML_ENKF_LOG_H
#ifdef __cplusplus
extern "C" {
#endif
#include <stdbool.h>
typedef struct rml_enkf_log_struct rml_enkf_log_type;
rml_enkf_log_type * rml_enkf_log_alloc();
void rml_enkf_log_free(rml_enkf_log_type * rml_log);
bool rml_enkf_log_get_clear_log( const rml_enkf_log_type * data );
void rml_enkf_log_set_clear_log( rml_enkf_log_type * data , bool clear_log);
void rml_enkf_log_set_log_file( rml_enkf_log_type * data , const char * log_file );
const char * rml_enkf_log_get_log_file( const rml_enkf_log_type * data);
void rml_enkf_log_open( rml_enkf_log_type * rml_log , int iteration_nr );
void rml_enkf_log_close( rml_enkf_log_type * rml_log );
void rml_enkf_log_line( rml_enkf_log_type * rml_log , const char * fmt , ...);
bool rml_enkf_log_is_open( const rml_enkf_log_type * rml_log );
#ifdef __cplusplus
}
#endif
#endif

148
ThirdParty/Ert/devel/libanalysis/src/analysis_module.c vendored
View File

@@ -1,19 +1,19 @@
/*
Copyright (C) 2011 Statoil ASA, Norway.
The file 'analysis_module.c' 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.
Copyright (C) 2011 Statoil ASA, Norway.
The file 'analysis_module.c' 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.
*/
#include <stdlib.h>
@@ -36,7 +36,7 @@ struct analysis_module_struct {
UTIL_TYPE_ID_DECLARATION;
void * lib_handle;
void * module_data;
char * symbol_table;
char * symbol_table;
char * lib_name;
analysis_free_ftype * freef;
@@ -49,17 +49,17 @@ struct analysis_module_struct {
analysis_get_options_ftype * get_options;
analysis_set_int_ftype * set_int;
analysis_set_double_ftype * set_double;
analysis_set_bool_ftype * set_bool;
analysis_set_bool_ftype * set_bool;
analysis_set_string_ftype * set_string;
analysis_has_var_ftype * has_var;
analysis_get_int_ftype * get_int;
analysis_get_double_ftype * get_double;
analysis_get_bool_ftype * get_bool;
analysis_get_ptr_ftype * get_ptr;
bool internal;
char * user_name; /* String used to identify this module for the user; not used in
bool internal;
char * user_name; /* String used to identify this module for the user; not used in
the linking process. */
};
@@ -99,15 +99,15 @@ static bool analysis_module_internal_check( analysis_module_type * module ) {
}
static analysis_module_type * analysis_module_alloc__( rng_type * rng ,
const analysis_table_type * table ,
const char * symbol_table ,
const char * lib_name ,
const char * user_name ,
static analysis_module_type * analysis_module_alloc__( rng_type * rng ,
const analysis_table_type * table ,
const char * symbol_table ,
const char * lib_name ,
const char * user_name ,
void * lib_handle ) {
analysis_module_type * module = analysis_module_alloc_empty( user_name , symbol_table , lib_name );
module->lib_handle = lib_handle;
module->initX = table->initX;
module->updateA = table->updateA;
@@ -119,7 +119,7 @@ static analysis_module_type * analysis_module_alloc__( rng_type * rng ,
module->set_bool = table->set_bool;
module->alloc = table->alloc;
module->freef = table->freef;
module->get_options = table->get_options;
module->get_options = table->get_options;
module->has_var = table->has_var;
module->get_int = table->get_int;
module->get_double = table->get_double;
@@ -142,11 +142,11 @@ static analysis_module_type * analysis_module_alloc__( rng_type * rng ,
static analysis_module_type * analysis_module_alloc( rng_type * rng ,
const char * user_name ,
const char * libname ,
static analysis_module_type * analysis_module_alloc( rng_type * rng ,
const char * user_name ,
const char * libname ,
const char * table_name ,
bool verbose,
bool verbose,
analysis_module_load_status_enum * load_status) {
analysis_module_type * module = NULL;
void * lib_handle = dlopen( libname , RTLD_NOW );
@@ -160,15 +160,15 @@ static analysis_module_type * analysis_module_alloc( rng_type * rng ,
if (verbose)
fprintf(stderr , "Failed to load symbol table:%s Error:%s \n",table_name , dlerror());
}
if (module == NULL)
if (module == NULL)
dlclose( lib_handle );
} else {
*load_status = DLOPEN_FAILURE;
if (verbose)
fprintf(stderr , "Failed to load library:%s Error:%s \n",libname , dlerror());
}
if (module != NULL) {
if (libname == NULL)
module->internal = true;
@@ -226,7 +226,7 @@ UTIL_IS_INSTANCE_FUNCTION( analysis_module , ANALYSIS_MODULE_TYPE_ID )
void analysis_module_free( analysis_module_type * module ) {
if (module->freef != NULL)
module->freef( module->module_data );
util_safe_free( module->lib_name );
free( module->user_name );
free( module->symbol_table );
@@ -244,26 +244,26 @@ void analysis_module_free__( void * arg) {
/* Update functions */
void analysis_module_initX(analysis_module_type * module ,
matrix_type * X ,
matrix_type * A ,
matrix_type * S ,
matrix_type * R ,
matrix_type * dObs ,
matrix_type * E ,
void analysis_module_initX(analysis_module_type * module ,
matrix_type * X ,
matrix_type * A ,
matrix_type * S ,
matrix_type * R ,
matrix_type * dObs ,
matrix_type * E ,
matrix_type * D ) {
module->initX(module->module_data , X , A , S , R , dObs , E , D );
}
void analysis_module_updateA(analysis_module_type * module ,
matrix_type * A ,
matrix_type * S ,
matrix_type * R ,
matrix_type * dObs ,
matrix_type * E ,
void analysis_module_updateA(analysis_module_type * module ,
matrix_type * A ,
matrix_type * S ,
matrix_type * R ,
matrix_type * dObs ,
matrix_type * E ,
matrix_type * D ) {
module->updateA(module->module_data , A , S , R , dObs , E , D );
@@ -272,12 +272,12 @@ void analysis_module_updateA(analysis_module_type * module ,
void analysis_module_init_update( analysis_module_type * module ,
const bool_vector_type * ens_mask ,
const matrix_type * S ,
const matrix_type * R ,
const matrix_type * dObs ,
const matrix_type * E ,
void analysis_module_init_update( analysis_module_type * module ,
const bool_vector_type * ens_mask ,
const matrix_type * S ,
const matrix_type * R ,
const matrix_type * dObs ,
const matrix_type * E ,
const matrix_type * D ) {
if (module->init_update != NULL)
module->init_update( module->module_data , ens_mask , S , R , dObs , E , D);
@@ -296,7 +296,7 @@ void analysis_module_complete_update( analysis_module_type * module ) {
static bool analysis_module_set_int(analysis_module_type * module , const char * flag , int value) {
if (module->set_int != NULL)
if (module->set_int != NULL)
return module->set_int( module->module_data , flag , value );
else
return false;
@@ -327,7 +327,7 @@ static bool analysis_module_set_string(analysis_module_type * module , const cha
/*
/*
The input value typically comes from the configuration system and
is in terms of a string, irrespective of the fundamental type of
the underlying parameter. The algorithm for setting the parameter
@@ -346,38 +346,38 @@ bool analysis_module_set_var( analysis_module_type * module , const char * var_n
bool set_ok = false;
{
int int_value;
if (util_sscanf_int( string_value , &int_value ))
if (util_sscanf_int( string_value , &int_value ))
set_ok = analysis_module_set_int( module , var_name , int_value );
if (set_ok)
return true;
}
{
double double_value;
if (util_sscanf_double( string_value , &double_value ))
if (util_sscanf_double( string_value , &double_value ))
set_ok = analysis_module_set_double( module , var_name , double_value );
if (set_ok)
return true;
}
{
bool bool_value;
if (util_sscanf_bool( string_value , &bool_value))
if (util_sscanf_bool( string_value , &bool_value))
set_ok = analysis_module_set_bool( module , var_name , bool_value );
if (set_ok)
return true;
}
set_ok = analysis_module_set_string( module , var_name , string_value );
if (!set_ok)
fprintf(stderr,"** Warning: failed to set %s=%s for analysis module:%s\n", var_name , string_value , module->user_name);
return set_ok;
}
@@ -404,7 +404,7 @@ int analysis_module_get_int( const analysis_module_type * module , const char *
return module->get_int( module->module_data , var );
else
util_exit("%s: Tried to get integer variable:%s from module:%s - get_int() method not implemented for this module\n" , __func__ , var , module->user_name);
} else
} else
util_exit("%s: Tried to get integer variable:%s from module:%s - module does not support this variable \n" , __func__ , var , module->user_name);
return 0;
@@ -417,7 +417,7 @@ bool analysis_module_get_bool( const analysis_module_type * module , const char
return module->get_bool( module->module_data , var );
else
util_exit("%s: Tried to get bool variable:%s from module:%s - get_int() method not implemented for this module\n" , __func__ , var , module->user_name);
} else
} else
util_exit("%s: Tried to get bool variable:%s from module:%s - module does not support this variable \n" , __func__ , var , module->user_name);
return false;
@@ -430,7 +430,7 @@ double analysis_module_get_double( const analysis_module_type * module , const c
return module->get_double( module->module_data , var );
else
util_exit("%s: Tried to get double variable:%s from module:%s - get_double() method not implemented for this module\n" , __func__ , var , module->user_name);
} else
} else
util_exit("%s: Tried to get double variable:%s from module:%s - module does not support this variable \n" , __func__ , var , module->user_name);
return 0;
@@ -443,9 +443,9 @@ void * analysis_module_get_ptr( const analysis_module_type * module , const char
return module->get_ptr( module->module_data , var );
else
util_exit("%s: Tried to get pointer variable:%s from module:%s - get_ptr() method not implemented for this module\n" , __func__ , var , module->user_name);
} else
} else
util_exit("%s: Tried to get pointer variable:%s from module:%s - module does not support this variable \n" , __func__ , var , module->user_name);
return NULL;
}

260
ThirdParty/Ert/devel/libanalysis/src/enkf_linalg.c vendored
View File

@@ -21,10 +21,10 @@ void enkf_linalg_genX3(matrix_type * X3 , const matrix_type * W , const matrix_t
for (i=0; i < nrmin; i++)
for (j=0; j < nrobs; j++)
matrix_iset(X1 , i , j , eig[i] * matrix_iget(W , j , i));
matrix_matmul(X2 , X1 , D); /* X2 = X1 * D (Eq. 14.31) */
matrix_matmul(X3 , W , X2); /* X3 = W * X2 = X1 * X2 (Eq. 14.31) */
matrix_matmul(X3 , W , X2); /* X3 = W * X2 = X1 * X2 (Eq. 14.31) */
matrix_free( X1 );
matrix_free( X2 );
}
@@ -34,7 +34,7 @@ void enkf_linalg_genX2(matrix_type * X2 , const matrix_type * S , const matrix_t
const int nrens = matrix_get_columns( S );
const int idim = matrix_get_rows( X2 );
matrix_dgemm(X2 , W , S , true , false , 1.0 , 0.0);
{
{
int i,j;
for (j=0; j < nrens; j++)
for (i=0; i < idim; i++)
@@ -53,11 +53,11 @@ void enkf_linalg_genX2(matrix_type * X2 , const matrix_type * S , const matrix_t
ncomp > 0 , truncation < 0: Use ncomp parameters.
ncomp < 0 , truncation > 0: Truncate at level 'truncation'.
The singular values are returned in the inv_sig0 vector; the values
we retain are inverted and the remaining elements in are explicitly
set to zero.
The left-hand singular vectors are returned in the matrix
U0. Depending on the value of the flag @store_V0T the right hand
singular vectors are stored in the V0T matrix, or just
@@ -74,14 +74,14 @@ void enkf_linalg_genX2(matrix_type * X2 , const matrix_type * S , const matrix_t
/*This function is similar to enkf_linalg_svdS but it returns the eigen values without its inverse and also give the matrices truncated U VT and Sig0*/
// Trunc.SVD(S) = U0 * Sig0 * V0T
int enkf_linalg_svd_truncation(const matrix_type * S ,
double truncation ,
int enkf_linalg_svd_truncation(const matrix_type * S ,
double truncation ,
int ncomp ,
dgesvd_vector_enum store_V0T ,
double * sig0,
matrix_type * U0 ,
dgesvd_vector_enum store_V0T ,
double * sig0,
matrix_type * U0 ,
matrix_type * V0T) {
int num_significant = -1;
int nrows = matrix_get_rows(S);
int ncolumns= matrix_get_columns(S);
@@ -90,9 +90,9 @@ int enkf_linalg_svd_truncation(const matrix_type * S ,
((truncation < 0) && (ncomp > 0))) {
int num_singular_values = util_int_min( matrix_get_rows( S ) , matrix_get_columns( S ));
{
{
matrix_type * workS = matrix_alloc_copy( S );
matrix_dgesvd(DGESVD_MIN_RETURN , store_V0T , workS , sig0 , U0 , V0T);
matrix_dgesvd(DGESVD_MIN_RETURN , store_V0T , workS , sig0 , U0 , V0T);
matrix_free( workS );
}
int i;
@@ -103,11 +103,11 @@ int enkf_linalg_svd_truncation(const matrix_type * S ,
double total_sigma2 = 0;
for (i=0; i < num_singular_values; i++)
total_sigma2 += sig0[i];
/*
/*
Determine the number of singular values by enforcing that
less than a fraction @truncation of the total variance be
accounted for.
accounted for.
*/
num_significant = 0;
{
@@ -116,40 +116,44 @@ int enkf_linalg_svd_truncation(const matrix_type * S ,
if (running_sigma2 / total_sigma2 < truncation) { /* Include one more singular value ? */
num_significant++;
running_sigma2 += sig0[i];
} else
} else
break;
}
}
}
matrix_resize(U0 , nrows , num_significant , true);
matrix_resize(V0T , num_significant , ncolumns , true);
if (num_significant > 0) {
matrix_resize(U0 , nrows , num_significant , true);
matrix_resize(V0T , num_significant , ncolumns , true);
} else
util_abort("%s: zero significant singular values\n",__func__);
}
else
else
util_abort("%s: truncation:%g ncomp:%d - invalid ambigous input.\n",__func__ , truncation , ncomp );
return num_significant;
}
int enkf_linalg_svdS(const matrix_type * S ,
double truncation ,
int enkf_linalg_svdS(const matrix_type * S ,
double truncation ,
int ncomp ,
dgesvd_vector_enum store_V0T ,
double * inv_sig0,
matrix_type * U0 ,
dgesvd_vector_enum store_V0T ,
double * inv_sig0,
matrix_type * U0 ,
matrix_type * V0T) {
double * sig0 = inv_sig0;
int num_significant = 0;
if (((truncation > 0) && (ncomp < 0)) ||
((truncation < 0) && (ncomp > 0))) {
int num_singular_values = util_int_min( matrix_get_rows( S ) , matrix_get_columns( S ));
{
{
matrix_type * workS = matrix_alloc_copy( S );
matrix_dgesvd(DGESVD_MIN_RETURN , store_V0T , workS , sig0 , U0 , V0T);
matrix_dgesvd(DGESVD_MIN_RETURN , store_V0T , workS , sig0 , U0 , V0T);
matrix_free( workS );
}
int i;
@@ -158,13 +162,13 @@ int enkf_linalg_svdS(const matrix_type * S ,
num_significant = ncomp;
else {
double total_sigma2 = 0;
for (i=0; i < num_singular_values; i++)
for (i=0; i < num_singular_values; i++)
total_sigma2 += sig0[i] * sig0[i];
/*
/*
Determine the number of singular values by enforcing that
less than a fraction @truncation of the total variance be
accounted for.
accounted for.
*/
num_significant = 0;
{
@@ -173,7 +177,7 @@ int enkf_linalg_svdS(const matrix_type * S ,
if (running_sigma2 / total_sigma2 < truncation) { /* Include one more singular value ? */
num_significant++;
running_sigma2 += sig0[i] * sig0[i];
} else
} else
break;
}
}
@@ -185,9 +189,9 @@ int enkf_linalg_svdS(const matrix_type * S ,
/* Explicitly setting the insignificant singular values to zero. */
for (i=num_significant; i < num_singular_values; i++)
inv_sig0[i] = 0;
} else
inv_sig0[i] = 0;
} else
util_abort("%s: truncation:%g ncomp:%d - invalid ambigous input.\n",__func__ , truncation , ncomp );
return num_significant;
}
@@ -200,14 +204,14 @@ void enkf_linalg_Cee(matrix_type * B, int nrens , const matrix_type * R , const
matrix_dgemm(X0 , U0 , R , true , false , 1.0 , 0.0); /* X0 = U0^T * R */
matrix_dgemm(B , X0 , U0 , false , false , 1.0 , 0.0); /* B = X0 * U0 */
matrix_free( X0 );
}
}
{
int i ,j;
/* Funny code ??
/* Funny code ??
Multiply B with S^(-1)from left and right
BHat = S^(-1) * B * S^(-1)
BHat = S^(-1) * B * S^(-1)
*/
for (j=0; j < matrix_get_columns( B ) ; j++)
for (i=0; i < matrix_get_rows( B ); i++)
@@ -217,26 +221,26 @@ void enkf_linalg_Cee(matrix_type * B, int nrens , const matrix_type * R , const
for (i=0; i < matrix_get_rows( B ); i++)
matrix_imul(B , i , j , inv_sig0[j]);
}
matrix_scale(B , nrens - 1.0);
}
void enkf_linalg_lowrankCinv__(const matrix_type * S ,
const matrix_type * R ,
matrix_type * V0T ,
matrix_type * Z,
double * eig ,
matrix_type * U0,
double truncation,
void enkf_linalg_lowrankCinv__(const matrix_type * S ,
const matrix_type * R ,
matrix_type * V0T ,
matrix_type * Z,
double * eig ,
matrix_type * U0,
double truncation,
int ncomp) {
const int nrobs = matrix_get_rows( S );
const int nrens = matrix_get_columns( S );
const int nrmin = util_int_min( nrobs , nrens );
double * inv_sig0 = util_calloc( nrmin , sizeof * inv_sig0);
if (V0T != NULL)
@@ -246,18 +250,18 @@ void enkf_linalg_lowrankCinv__(const matrix_type * S ,
{
matrix_type * B = matrix_alloc( nrmin , nrmin );
enkf_linalg_Cee( B , nrens , R , U0 , inv_sig0); /* B = Xo = (N-1) * Sigma0^(+) * U0'* Cee * U0 * Sigma0^(+') (14.26)*/
enkf_linalg_Cee( B , nrens , R , U0 , inv_sig0); /* B = Xo = (N-1) * Sigma0^(+) * U0'* Cee * U0 * Sigma0^(+') (14.26)*/
matrix_dgesvd(DGESVD_MIN_RETURN , DGESVD_NONE, B , eig, Z , NULL);
matrix_free( B );
}
{
int i,j;
/* Lambda1 = (I + Lambda)^(-1) */
for (i=0; i < nrmin; i++)
for (i=0; i < nrmin; i++)
eig[i] = 1.0 / (1 + eig[i]);
for (j=0; j < nrmin; j++)
for (i=0; i < nrmin; i++)
matrix_imul(Z , i , j , inv_sig0[i]); /* Z2 = Sigma0^(+) * Z; */
@@ -266,20 +270,20 @@ void enkf_linalg_lowrankCinv__(const matrix_type * S ,
}
void enkf_linalg_lowrankCinv(const matrix_type * S ,
const matrix_type * R ,
void enkf_linalg_lowrankCinv(const matrix_type * S ,
const matrix_type * R ,
matrix_type * W , /* Corresponding to X1 from Eq. 14.29 */
double * eig , /* Corresponding to 1 / (1 + Lambda_1) (14.29) */
double truncation ,
int ncomp) {
const int nrobs = matrix_get_rows( S );
const int nrens = matrix_get_columns( S );
const int nrmin = util_int_min( nrobs , nrens );
matrix_type * U0 = matrix_alloc( nrobs , nrmin );
matrix_type * Z = matrix_alloc( nrmin , nrmin );
enkf_linalg_lowrankCinv__( S , R , NULL , Z , eig , U0 , truncation , ncomp);
matrix_matmul(W , U0 , Z); /* X1 = W = U0 * Z2 = U0 * Sigma0^(+') * Z */
@@ -288,13 +292,13 @@ void enkf_linalg_lowrankCinv(const matrix_type * S ,
}
void enkf_linalg_meanX5(const matrix_type * S ,
const matrix_type * W ,
const double * eig ,
const matrix_type * dObs,
void enkf_linalg_meanX5(const matrix_type * S ,
const matrix_type * W ,
const double * eig ,
const matrix_type * dObs,
matrix_type * X5) {
const int nrens = matrix_get_columns( S );
const int nrobs = matrix_get_rows( S );
const int nrmin = util_int_min( nrobs , nrens );
@@ -304,8 +308,8 @@ void enkf_linalg_meanX5(const matrix_type * S ,
double * y1 = &work[0];
double * y2 = &work[nrmin];
double * y3 = &work[2*nrmin];
double * y4 = &work[2*nrmin + nrobs];
double * y4 = &work[2*nrmin + nrobs];
if (nrobs == 1) {
/* Is this special casing necessary ??? */
y1[0] = matrix_iget(W , 0,0) * matrix_iget( innov , 0 , 0);
@@ -317,13 +321,13 @@ void enkf_linalg_meanX5(const matrix_type * S ,
matrix_dgemv(W , matrix_get_data( innov ) , y1 , true , 1.0, 0.0); /* y1 = Trans(W) * innov */
for (int i= 0; i < nrmin; i++)
y2[i] = eig[i] * y1[i]; /* y2 = eig * y1 */
matrix_dgemv(W , y2 , y3 , false , 1.0 , 0.0); /* y3 = W * y2; */
matrix_dgemv(W , y2 , y3 , false , 1.0 , 0.0); /* y3 = W * y2; */
matrix_dgemv(S , y3 , y4 , true , 1.0 , 0.0); /* y4 = Trans(S) * y3 */
}
for (int iens = 0; iens < nrens; iens++)
matrix_set_column(X5 , y4 , iens );
matrix_shift(X5 , 1.0/nrens);
}
free( work );
@@ -332,7 +336,7 @@ void enkf_linalg_meanX5(const matrix_type * S ,
void enkf_linalg_X5sqrt(matrix_type * X2 , matrix_type * X5 , const matrix_type * randrot, int nrobs) {
void enkf_linalg_X5sqrt(matrix_type * X2 , matrix_type * X5 , const matrix_type * randrot, int nrobs) {
const int nrens = matrix_get_columns( X5 );
const int nrmin = util_int_min( nrobs , nrens );
matrix_type * VT = matrix_alloc( nrens , nrens );
@@ -348,24 +352,24 @@ void enkf_linalg_X5sqrt(matrix_type * X2 , matrix_type * X5 , const matrix_type
int i,j;
for (i = 0; i < nrmin; i++)
isig[i] = sqrt( util_double_max( 1.0 - sig[i]*sig[i] ,0.0));
for (j = 0; j < nrens; j++)
for (i = 0; i < nrens; i++)
matrix_iset(X3 , i , j , matrix_iget(VT , j , i));
for (j=0; j< nrmin; j++)
matrix_scale_column(X3 , j , isig[j]);
matrix_dgemm(X33 , X3 , VT , false , false , 1.0 , 0.0); /* X33 = X3 * VT */
if (randrot != NULL)
matrix_dgemm(X4 , X33 , randrot , false, false , 1.0 , 0.0); /* X4 = X33 * Randrot */
matrix_dgemm(X4 , X33 , randrot , false, false , 1.0 , 0.0); /* X4 = X33 * Randrot */
else
matrix_assign(X4 , X33);
matrix_set(IenN , -1.0/ nrens);
for (i = 0; i < nrens; i++)
matrix_iadd(IenN , i , i , 1.0);
matrix_dgemm(X5 , IenN , X4 , false , false , 1.0 , 1.0); /* X5 = IenN * X4 + X5 */
matrix_free( X3 );
@@ -383,51 +387,51 @@ void enkf_linalg_X5sqrt(matrix_type * X2 , matrix_type * X5 , const matrix_type
matrix_type * enkf_linalg_alloc_innov( const matrix_type * dObs , const matrix_type * S) {
matrix_type * innov = matrix_alloc_copy( dObs );
for (int iobs =0; iobs < matrix_get_row_sum( dObs , iobs); iobs++)
for (int iobs =0; iobs < matrix_get_row_sum( dObs , iobs); iobs++)
matrix_isub( innov , iobs , 0 , matrix_get_row_sum( S , iobs ));
return innov;
}
void enkf_linalg_init_stdX( matrix_type * X , const matrix_type * S , const matrix_type * D ,
void enkf_linalg_init_stdX( matrix_type * X , const matrix_type * S , const matrix_type * D ,
const matrix_type * W , const double * eig , bool bootstrap) {
int nrobs = matrix_get_rows( W );
int ens_size = matrix_get_rows( X );
matrix_type * X3 = matrix_alloc(nrobs , ens_size);
enkf_linalg_genX3(X3 , W , D , eig ); /* X2 = diag(eig) * W' * D (Eq. 14.31, Evensen (2007)) */
/* X3 = W * X2 = X1 * X2 (Eq. 14.31, Evensen (2007)) */
/* X3 = W * X2 = X1 * X2 (Eq. 14.31, Evensen (2007)) */
matrix_dgemm( X , S , X3 , true , false , 1.0 , 0.0); /* X = S' * X3 */
if (!bootstrap) {
for (int i = 0; i < ens_size ; i++)
matrix_iadd( X , i , i , 1.0); /*X = I + X */
}
matrix_free( X3 );
}
void enkf_linalg_init_sqrtX(matrix_type * X5 ,
const matrix_type * S ,
const matrix_type * randrot ,
const matrix_type * innov ,
const matrix_type * W ,
const double * eig ,
void enkf_linalg_init_sqrtX(matrix_type * X5 ,
const matrix_type * S ,
const matrix_type * randrot ,
const matrix_type * innov ,
const matrix_type * W ,
const double * eig ,
bool bootstrap) {
const int nrobs = matrix_get_rows( S );
const int nrens = matrix_get_columns( S );
const int nrmin = util_int_min( nrobs , nrens );
matrix_type * X2 = matrix_alloc(nrmin , nrens);
if (bootstrap)
util_exit("%s: Sorry bootstrap support not fully implemented for SQRT scheme\n",__func__);
enkf_linalg_meanX5( S , W , eig , innov , X5 );
enkf_linalg_genX2(X2 , S , W , eig);
enkf_linalg_X5sqrt(X2 , X5 , randrot , nrobs);
@@ -450,14 +454,14 @@ void enkf_linalg_init_sqrtX(matrix_type * X5 ,
2. Then the QR decomposition is computed, and Q will then be a random orthogonal matrix.
3. The diagonal elements of R are extracted and we construct the diagonal matrix X(j,j)=R(j,j)/|R(j,j)|
4. An updated Q'=Q X is computed, and this is now a random orthogonal matrix with a Haar measure.
The implementation is a plain reimplementation/copy of the old m_randrot.f90 function.
*/
/**
NB: This should rather use the implementation in m_mean_preserving_rotation.f90.
NB: This should rather use the implementation in m_mean_preserving_rotation.f90.
*/
void enkf_linalg_set_randrot( matrix_type * Q , rng_type * rng) {
@@ -465,8 +469,8 @@ void enkf_linalg_set_randrot( matrix_type * Q , rng_type * rng) {
double * tau = util_calloc( ens_size , sizeof * tau );
int * sign = util_calloc( ens_size , sizeof * sign);
for (int i = 0; i < ens_size; i++)
for (int j = 0; j < ens_size; j++)
for (int i = 0; i < ens_size; i++)
for (int j = 0; j < ens_size; j++)
matrix_iset(Q , i , j , rng_std_normal( rng ));
matrix_dgeqrf( Q , tau ); /* QR factorization */
@@ -482,7 +486,7 @@ void enkf_linalg_set_randrot( matrix_type * Q , rng_type * rng) {
}
free( sign );
free( tau );
free( tau );
}
@@ -492,7 +496,7 @@ void enkf_linalg_set_randrot( matrix_type * Q , rng_type * rng) {
Generates the mean preserving random rotation for the EnKF SQRT algorithm
using the algorithm from Sakov 2006-07. I.e, generate rotation Up such that
Up*Up^T=I and Up*1=1 (all rows have sum = 1) see eq 17.
From eq 18, Up=B * Upb * B^T
From eq 18, Up=B * Upb * B^T
B is a random orthonormal basis with the elements in the first column equals 1/sqrt(nrens)
Upb = | 1 0 |
@@ -507,8 +511,8 @@ matrix_type * enkf_linalg_alloc_mp_randrot(int ens_size , rng_type * rng) {
matrix_type * B = matrix_alloc( ens_size , ens_size );
matrix_type * Upb = matrix_alloc( ens_size , ens_size );
matrix_type * U = matrix_alloc_shared(Upb , 1 , 1 , ens_size - 1, ens_size - 1);
{
int k,j;
matrix_type * R = matrix_alloc( ens_size , ens_size );
@@ -517,7 +521,7 @@ matrix_type * enkf_linalg_alloc_mp_randrot(int ens_size , rng_type * rng) {
/* modified_gram_schmidt is used to create the orthonormal basis in B.*/
for (k=0; k < ens_size; k++) {
double Rkk = sqrt( matrix_column_column_dot_product( B , k , B , k));
double Rkk = sqrt( matrix_column_column_dot_product( B , k , B , k));
matrix_iset(R , k , k , Rkk);
matrix_scale_column(B , k , 1.0/Rkk);
for (j=k+1; j < ens_size; j++) {
@@ -535,23 +539,23 @@ matrix_type * enkf_linalg_alloc_mp_randrot(int ens_size , rng_type * rng) {
}
matrix_free( R );
}
enkf_linalg_set_randrot( U , rng );
matrix_iset( Upb , 0 , 0 , 1);
{
matrix_type * Q = matrix_alloc( ens_size , ens_size );
matrix_dgemm( Q , B , Upb , false , false , 1, 0); /* Q = B * Ubp */
matrix_dgemm( Up , Q , B , false , true , 1, 0); /* Up = Q * T(B) */
matrix_free( Q );
}
matrix_free( B );
matrix_free( Upb );
matrix_free( U );
}
return Up;
}
@@ -561,7 +565,7 @@ matrix_type * enkf_linalg_alloc_mp_randrot(int ens_size , rng_type * rng) {
/**
Checking that the sum through one row in the X matrix equals
@target_sum. @target_sum will be 1 normally, and zero if we are doing
bootstrap.
bootstrap.
*/
void enkf_linalg_checkX(const matrix_type * X , bool bootstrap) {
@@ -572,10 +576,10 @@ void enkf_linalg_checkX(const matrix_type * X , bool bootstrap) {
target_sum = 0;
else
target_sum = 1;
for (int icol = 0; icol < matrix_get_columns( X ); icol++) {
double col_sum = matrix_get_column_sum(X , icol);
if (fabs(col_sum - target_sum) > 0.0001)
if (fabs(col_sum - target_sum) > 0.0001)
util_abort("%s: something is seriously broken. col:%d col_sum = %g != %g - ABORTING\n",__func__ , icol , col_sum , target_sum);
}
}
@@ -584,14 +588,14 @@ void enkf_linalg_checkX(const matrix_type * X , bool bootstrap) {
/*****************************************************************/
void enkf_linalg_get_PC( const matrix_type * S0,
const matrix_type * dObs ,
void enkf_linalg_get_PC( const matrix_type * S0,
const matrix_type * dObs ,
double truncation,
int ncomp,
int ncomp,
matrix_type * PC,
matrix_type * PC_obs,
matrix_type * PC_obs,
double_vector_type * singular_values) {
const int nrobs = matrix_get_rows( S0 );
const int nrens = matrix_get_columns( S0 );
const int nrmin = util_int_min( nrobs , nrens );
@@ -607,7 +611,7 @@ void enkf_linalg_get_PC( const matrix_type * S0,
{
matrix_type * S_mean = matrix_alloc( nrobs , 1 );
int num_PC = enkf_linalg_svdS(S , truncation , ncomp, DGESVD_NONE , inv_sig0 , U0 , NULL);
matrix_assign( S , S0); // The svd routine will overwrite S - we therefor must pick it up again from S0.
matrix_subtract_and_store_row_mean( S , S_mean);
@@ -622,7 +626,7 @@ void enkf_linalg_get_PC( const matrix_type * S0,
matrix_dgemm( PC , U0 , S , true , false , 1.0 , 0.0 );
}
/* The observer projections. */
{
matrix_scale( S_mean , -1.0);
@@ -633,7 +637,7 @@ void enkf_linalg_get_PC( const matrix_type * S0,
for (int i=0; i < double_vector_size( singular_values ); i++)
inv_sig0[i] = 1.0 / inv_sig0[i];
matrix_free( S_mean );
}
matrix_free( S );
@@ -649,7 +653,7 @@ void enkf_linalg_rml_enkfX1(matrix_type *X1, matrix_type *Udr, matrix_type *D, m
here the negative sign cancels with one needed in X3 matrix computation
*/
matrix_type * tmp = matrix_alloc(matrix_get_columns(Udr),matrix_get_rows(R));
matrix_matmul_with_transpose( tmp, Udr, R, true, false);
matrix_matmul( X1 , tmp, D);
@@ -661,23 +665,23 @@ void enkf_linalg_rml_enkfX1(matrix_type *X1, matrix_type *Udr, matrix_type *D, m
void enkf_linalg_rml_enkfX2(matrix_type *X2 , double *Wdr , matrix_type * X1 , double a , int nsign)
{
/*
This routine computes X2 for RML_EnKF module as X2 = ((a*Ipd)+Wd^2)^-1 * X1
This routine computes X2 for RML_EnKF module as X2 = ((a*Ipd)+Wd^2)^-1 * X1
Since a+Ipd & Wd are diagonal in nature the computation is reduced to array operations
*/
for (int i=0; i< nsign ; i++) {
double scale_factor = 1 / (a + (Wdr[i]*Wdr[i]));
matrix_scale_row(X1 , i , scale_factor);
}
matrix_assign(X2,X1);
}
void enkf_linalg_rml_enkfX3(matrix_type *X3, matrix_type *VdTr, double *Wdr, matrix_type *X2, int nsign)
{
/*
This routine computes X3 for RML_EnKF module as X3 = Vd *Wd*X2
This routine computes X3 for RML_EnKF module as X3 = Vd *Wd*X2
*/
printf("\nWd: ");
@@ -701,9 +705,9 @@ double enkf_linalg_data_mismatch(matrix_type *D , matrix_type *R , matrix_type *
double mismatch;
matrix_matmul_with_transpose(tmp, D, R,true, false); // tmp = D' * R, i.e. N-by-p
matrix_matmul(Sk, tmp, D); // Sk = D' * R * D
// Calculate the mismatch
// Calculate the mismatch
mismatch = matrix_trace(Sk)/(matrix_get_columns(D));
return mismatch;
}
@@ -718,7 +722,7 @@ void enkf_linalg_Covariance(matrix_type *Cd, const matrix_type *E, double nsc ,i
}
}
nsc = nsc*nsc;
matrix_scale(Cd,nsc);
matrix_scale(Cd,nsc);
}
@@ -746,7 +750,7 @@ void enkf_linalg_rml_enkfX7(matrix_type * X7, matrix_type * VdT, double * Wdr, d
double scale_factor = 1 / ( a + (Wdr[i]*Wdr[i]));
matrix_scale_row( tmp1 , i , scale_factor);
}
matrix_matmul_with_transpose(tmp2, tmp1, VdT, true, false);
matrix_matmul(X7, tmp2, X6);