Commit dde3f379 authored by sebastien's avatar sebastien
Browse files

Build system:

* Simplified the way we deal with various Octave/MATLAB contexts from MEX files:
  - only three defines: MATLAB_MEX_FILE, MATLAB_VERSION (hex number) and OCTAVE_MEX_FILE
  - one header for MEX files: dynmex.h
  - headers for BLAS and LAPACK: dynblas.h and dynlapack.h (used from Dynare++ and the MEX files)
* Merged the two sources trees of sylvester library


git-svn-id: https://www.dynare.org/svn/dynare/trunk@3006 ac1d8469-bf42-47a9-8791-bf33cf982152
parent 834c6002
......@@ -25,7 +25,7 @@ GENERATED_FILES = \
noinst_LIBRARIES = libinteg.a
libinteg_a_SOURCES = $(CWEBSRC) $(GENERATED_FILES) precalc_quadrature.dat
libinteg_a_CPPFLAGS = -I../../sylv/cc -I../../tl/cc -DPOSIX_THREADS
libinteg_a_CPPFLAGS = -I../../sylv/cc -I../../tl/cc -I$(top_srcdir)/mex/sources -DPOSIX_THREADS
libinteg_a_CXXFLAGS = $(PTHREAD_CFLAGS)
BUILT_SOURCES = $(GENERATED_FILES)
......
......@@ -6,7 +6,8 @@
@c
#include "vector_function.h"
#include "cpplapack.h"
#include <dynlapack.h>
#include <math.h>
......@@ -169,7 +170,7 @@ void GaussConverterFunction::calcCholeskyFactor(const GeneralMatrix& vcov)
A.get(i,j) = 0.0;
int info;
LAPACK_dpotrf("L", &rows, A.base(), &rows, &info);
dpotrf("L", &rows, A.base(), &rows, &info);
// todo: raise if |info!=1|
}
......
check_PROGRAMS = tests
tests_SOURCES = tests.cpp
tests_CPPFLAGS = -I../cc -I../../tl/cc -I../../sylv/cc -DPOSIX_THREADS
tests_CPPFLAGS = -I../cc -I../../tl/cc -I../../sylv/cc -I$(top_srcdir)/mex/sources -DPOSIX_THREADS
tests_CXXFLAGS = $(PTHREAD_CFLAGS)
tests_LDADD = ../../tl/cc/libtl.a ../../sylv/cc/libsylv.a ../cc/libinteg.a $(LAPACK_LIBS) $(BLAS_LIBS) $(LIBS) $(FLIBS) $(PTHREAD_LIBS)
......
......@@ -2,7 +2,7 @@
/* Copyright 2005, Ondra Kamenik */
#include "GeneralMatrix.h"
#include "cpplapack.h"
#include <dynlapack.h>
#include "SylvException.h"
#include "rfs_tensor.h"
......@@ -260,7 +260,7 @@ bool TestRunnable::qmc_normal_moments(const GeneralMatrix& m, int imom, int leve
for (int j = i+1; j < rows; j++)
mchol.get(i,j) = 0.0;
int info;
LAPACK_dpotrf("L", &rows, mchol.base(), &rows, &info);
dpotrf("L", &rows, mchol.base(), &rows, &info);
// make vector function
MomentFunction func(mchol, imom);
......
......@@ -53,7 +53,7 @@ GENERATED_FILES = \
noinst_LIBRARIES = libkord.a
libkord_a_SOURCES = $(CWEBSRC) $(GENERATED_FILES)
libkord_a_CPPFLAGS = -I../sylv/cc -I../tl/cc -I../integ/cc -DPOSIX_THREADS
libkord_a_CPPFLAGS = -I../sylv/cc -I../tl/cc -I../integ/cc -I$(top_srcdir)/mex/sources -DPOSIX_THREADS
libkord_a_CXXFLAGS = $(PTHREAD_CFLAGS)
BUILT_SOURCES = $(GENERATED_FILES)
......
......@@ -9,7 +9,8 @@
#include "dynamic_model.h"
#include "SymSchurDecomp.h"
#include "cpplapack.h"
#include <dynlapack.h>
#include <limits>
......@@ -609,7 +610,7 @@ void RandomShockRealization::choleskyFactor(const TwoDMatrix& v)
factor.get(i,j) = 0.0;
int info;
LAPACK_dpotrf("L", &rows, factor.base(), &rows, &info);
dpotrf("L", &rows, factor.base(), &rows, &info);
KORD_RAISE_IF(info != 0,
"Info!=0 in RandomShockRealization::choleskyFactor");
}
......
......@@ -7,7 +7,8 @@
#include "kord_exception.h"
#include "first_order.h"
#include "cpplapack.h"
#include <dynlapack.h>
double qz_criterium = 1.000001;
@<|order_eigs| function code@>;
......@@ -15,7 +16,7 @@ double qz_criterium = 1.000001;
@<|FirstOrder::journalEigs| code@>;
@ This is a function which selects the eigenvalues pair used by
|LAPACK_dgges|. See documentation to DGGES for details. Here we want
|dgges|. See documentation to DGGES for details. Here we want
to select (return true) the pairs for which $\alpha<\beta$.
@<|order_eigs| function code@>=
......@@ -192,7 +193,7 @@ the same. The difference is only numerical error.
Vector work(lwork);
IntSequence bwork(n);
int info;
LAPACK_dgges("N", "V", "S", order_eigs, &n, matE.getData().base(), &n,
dgges("N", "V", "S", order_eigs, &n, matE.getData().base(), &n,
matD.getData().base(), &n, &sdim, alphar.base(), alphai.base(),
beta.base(), vsl.getData().base(), &n, vsr.getData().base(), &n,
work.base(), &lwork, &(bwork[0]), &info);
......
......@@ -12,8 +12,6 @@
#include "product.h"
#include "quasi_mcarlo.h"
#include "cpplapack.h"
@<|ResidFunction| constructor code@>;
@<|ResidFunction| copy constructor code@>;
@<|ResidFunction| destructor code@>;
......
......@@ -8,7 +8,7 @@
#include "kord_exception.h"
#include "korder.h"
#include "cpplapack.h"
#include <dynlapack.h>
@<|PLUMatrix| copy constructor@>;
@<|PLUMatrix::calcPLU| code@>;
......@@ -40,7 +40,7 @@ void PLUMatrix::calcPLU()
int info;
int rows = nrows();
inv = (const Vector&)getData();
LAPACK_dgetrf(&rows, &rows, inv.base(), &rows, ipiv, &info);
dgetrf(&rows, &rows, inv.base(), &rows, ipiv, &info);
}
@ Here we just call the LAPACK machinery to multiply by the inverse.
......@@ -54,7 +54,7 @@ void PLUMatrix::multInv(TwoDMatrix& m) const
int mcols = m.ncols();
int mrows = m.nrows();
double* mbase = m.getData().base();
LAPACK_dgetrs("N", &mrows, &mcols, inv.base(), &mrows, ipiv,
dgetrs("N", &mrows, &mcols, inv.base(), &mrows, ipiv,
mbase, &mrows, &info);
KORD_RAISE_IF(info != 0,
"Info!=0 in PLUMatrix::multInv");
......
......@@ -6,8 +6,8 @@
#include "SylvException.h"
#include "GeneralMatrix.h"
#include "cppblas.h"
#include "cpplapack.h"
#include <dynblas.h>
#include <dynlapack.h>
#include <stdio.h>
#include <string.h>
......@@ -270,14 +270,14 @@ void GeneralMatrix::gemm(const char* transa, const ConstGeneralMatrix& a,
throw SYLV_MES_EXCEPTION("Wrong dimensions for matrix multiplication.");
}
int m = opa_rows;
int n = opb_cols;
int k = opa_cols;
int lda = a.ld;
int ldb = b.ld;
int ldc = ld;
blas_int m = opa_rows;
blas_int n = opb_cols;
blas_int k = opa_cols;
blas_int lda = a.ld;
blas_int ldb = b.ld;
blas_int ldc = ld;
if (lda > 0 && ldb > 0 && ldc > 0) {
BLAS_dgemm(transa, transb, &m, &n, &k, &alpha, a.data.base(), &lda,
dgemm(transa, transb, &m, &n, &k, &alpha, a.data.base(), &lda,
b.data.base(), &ldb, &beta, data.base(), &ldc);
} else if (numRows()*numCols() > 0) {
if (beta == 0.0)
......@@ -365,14 +365,14 @@ void ConstGeneralMatrix::multVec(double a, Vector& x, double b, const ConstVecto
throw SYLV_MES_EXCEPTION("Wrong dimensions for vector multiply.");
}
if (rows > 0) {
int mm = rows;
int nn = cols;
blas_int mm = rows;
blas_int nn = cols;
double alpha = b;
int lda = ld;
int incx = d.skip();
blas_int lda = ld;
blas_int incx = d.skip();
double beta = a;
int incy = x.skip();
BLAS_dgemv("N", &mm, &nn, &alpha, data.base(), &lda, d.base(), &incx,
blas_int incy = x.skip();
dgemv("N", &mm, &nn, &alpha, data.base(), &lda, d.base(), &incx,
&beta, x.base(), &incy);
}
......@@ -385,14 +385,14 @@ void ConstGeneralMatrix::multVecTrans(double a, Vector& x, double b,
throw SYLV_MES_EXCEPTION("Wrong dimensions for vector multiply.");
}
if (rows > 0) {
int mm = rows;
int nn = cols;
blas_int mm = rows;
blas_int nn = cols;
double alpha = b;
int lda = rows;
int incx = d.skip();
blas_int lda = rows;
blas_int incx = d.skip();
double beta = a;
int incy = x.skip();
BLAS_dgemv("T", &mm, &nn, &alpha, data.base(), &lda, d.base(), &incx,
blas_int incy = x.skip();
dgemv("T", &mm, &nn, &alpha, data.base(), &lda, d.base(), &incx,
&beta, x.base(), &incy);
}
}
......@@ -409,11 +409,12 @@ void ConstGeneralMatrix::multInvLeft(const char* trans, int mrows, int mcols, in
if (rows > 0) {
GeneralMatrix inv(*this);
int* ipiv = new int[rows];
int info;
LAPACK_dgetrf(&rows, &rows, inv.getData().base(), &rows, ipiv, &info);
LAPACK_dgetrs(trans, &rows, &mcols, inv.base(), &rows, ipiv, d,
&mld, &info);
lapack_int* ipiv = new lapack_int[rows];
lapack_int info;
lapack_int rows2 = rows, mrows2 = mrows, mcols2 = mcols, mld2 = mld;
dgetrf(&rows2, &rows2, inv.getData().base(), &rows2, ipiv, &info);
dgetrs(trans, &rows2, &mcols2, inv.base(), &rows2, ipiv, d,
&mld2, &info);
delete [] ipiv;
}
}
......
noinst_LIBRARIES = libsylv.a
# For dynblas.h and dynlapack.h
libsylv_a_CPPFLAGS = -I$(top_srcdir)/mex/sources
libsylv_a_SOURCES = \
IterativeSylvester.cpp \
cpplapack.h \
SylvMatrix.h \
QuasiTriangular.cpp \
QuasiTriangularZero.cpp \
......@@ -19,7 +21,6 @@ libsylv_a_SOURCES = \
BlockDiagonal.h \
SylvesterSolver.h \
SylvException.cpp \
cppblas.h \
SimilarityDecomp.h \
IterativeSylvester.h \
SchurDecompEig.cpp \
......
......@@ -6,7 +6,7 @@
#include "SylvException.h"
#include "SchurDecomp.h"
#include "cppblas.h"
#include <dynblas.h>
#include <stdio.h>
#include <cmath>
......@@ -396,10 +396,10 @@ QuasiTriangular::QuasiTriangular(int p, const QuasiTriangular& t)
: SqSylvMatrix(t.numRows()), diagonal(getData().base(), t.diagonal)
{
Vector aux(t.getData());
int d_size = diagonal.getSize();
blas_int d_size = diagonal.getSize();
double alpha = 1.0;
double beta = 0.0;
BLAS_dgemm("N", "N", &d_size, &d_size, &d_size, &alpha, aux.base(),
dgemm("N", "N", &d_size, &d_size, &d_size, &alpha, aux.base(),
&d_size, t.getData().base(), &d_size, &beta, getData().base(), &d_size);
}
......@@ -529,10 +529,10 @@ void QuasiTriangular::solvePre(Vector& x, double& eig_min)
eig_min = eig_size;
}
int nn = diagonal.getSize();
int lda = diagonal.getSize();
int incx = x.skip();
BLAS_dtrsv("U", "N", "N", &nn, getData().base(), &lda, x.base(), &incx);
blas_int nn = diagonal.getSize();
blas_int lda = diagonal.getSize();
blas_int incx = x.skip();
dtrsv("U", "N", "N", &nn, getData().base(), &lda, x.base(), &incx);
}
void QuasiTriangular::solvePreTrans(Vector& x, double& eig_min)
......@@ -550,10 +550,10 @@ void QuasiTriangular::solvePreTrans(Vector& x, double& eig_min)
eig_min = eig_size;
}
int nn = diagonal.getSize();
int lda = diagonal.getSize();
int incx = x.skip();
BLAS_dtrsv("U", "T", "N", &nn, getData().base(), &lda, x.base(), &incx);
blas_int nn = diagonal.getSize();
blas_int lda = diagonal.getSize();
blas_int incx = x.skip();
dtrsv("U", "T", "N", &nn, getData().base(), &lda, x.base(), &incx);
}
......@@ -561,10 +561,10 @@ void QuasiTriangular::solvePreTrans(Vector& x, double& eig_min)
void QuasiTriangular::multVec(Vector& x, const ConstVector& b) const
{
x = b;
int nn = diagonal.getSize();
int lda = diagonal.getSize();
int incx = x.skip();
BLAS_dtrmv("U", "N", "N", &nn, getData().base(), &lda, x.base(), &incx);
blas_int nn = diagonal.getSize();
blas_int lda = diagonal.getSize();
blas_int incx = x.skip();
dtrmv("U", "N", "N", &nn, getData().base(), &lda, x.base(), &incx);
for (const_diag_iter di = diag_begin(); di != diag_end(); ++di) {
if (!(*di).isReal()) {
int jbar = (*di).getIndex();
......@@ -577,10 +577,10 @@ void QuasiTriangular::multVec(Vector& x, const ConstVector& b) const
void QuasiTriangular::multVecTrans(Vector& x, const ConstVector& b) const
{
x = b;
int nn = diagonal.getSize();
int lda = diagonal.getSize();
int incx = x.skip();
BLAS_dtrmv("U", "T", "N", &nn, getData().base(), &lda, x.base(), &incx);
blas_int nn = diagonal.getSize();
blas_int lda = diagonal.getSize();
blas_int incx = x.skip();
dtrmv("U", "T", "N", &nn, getData().base(), &lda, x.base(), &incx);
for (const_diag_iter di = diag_begin(); di != diag_end(); ++di) {
if (!(*di).isReal()) {
int jbar = (*di).getIndex();
......
......@@ -4,21 +4,21 @@
#include "SchurDecomp.h"
#include "cpplapack.h"
#include <dynlapack.h>
SchurDecomp::SchurDecomp(const SqSylvMatrix& m)
: q_destroy(true), t_destroy(true)
{
int rows = m.numRows();
lapack_int rows = m.numRows();
q = new SqSylvMatrix(rows);
SqSylvMatrix auxt(m);
int sdim;
lapack_int sdim;
double* const wr = new double[rows];
double* const wi = new double[rows];
int lwork = 6*rows;
lapack_int lwork = 6*rows;
double* const work = new double[lwork];
int info;
LAPACK_dgees("V", "N", 0, &rows, auxt.base(), &rows, &sdim,
lapack_int info;
dgees("V", "N", 0, &rows, auxt.base(), &rows, &sdim,
wr, wi, q->base(), &rows,
work, &lwork, 0, &info);
delete [] work;
......
#include "SchurDecompEig.h"
#include "SylvException.h"
#include "cpplapack.h"
#include <dynlapack.h>
/* bubble diagonal 1-1, or 2-2 block from position 'from' to position
* 'to'. If an eigenvalue cannot be swapped with its neighbour, the
......@@ -45,16 +46,16 @@ bool SchurDecompEig::tryToSwap(diag_iter& it, diag_iter& itadd)
itadd = it;
--itadd;
int n = getDim();
int ifst = (*it).getIndex() + 1;
int ilst = (*itadd).getIndex() + 1;
lapack_int n = getDim();
lapack_int ifst = (*it).getIndex() + 1;
lapack_int ilst = (*itadd).getIndex() + 1;
double* work = new double[n];
int info;
LAPACK_dtrexc("V", &n, getT().base(), &n, getQ().base(), &n, &ifst, &ilst, work,
lapack_int info;
dtrexc("V", &n, getT().base(), &n, getQ().base(), &n, &ifst, &ilst, work,
&info);
delete [] work;
if (info < 0) {
throw SYLV_MES_EXCEPTION("Wrong argument to LAPACK_dtrexc.");
throw SYLV_MES_EXCEPTION("Wrong argument to dtrexc.");
}
if (info == 0) {
......
......@@ -7,7 +7,7 @@
#include "SchurDecompEig.h"
#include "SylvException.h"
#include "cpplapack.h"
#include <dynlapack.h>
#include <cmath>
......@@ -53,12 +53,12 @@ bool SimilarityDecomp::solveX(diag_iter start, diag_iter end,
SqSylvMatrix B((const GeneralMatrix&)*b, ei, ei, X.numCols());
GeneralMatrix C((const GeneralMatrix&)*b, si, ei, X.numRows(), X.numCols());
int isgn = -1;
int m = A.numRows();
int n = B.numRows();
lapack_int isgn = -1;
lapack_int m = A.numRows();
lapack_int n = B.numRows();
double scale;
int info;
LAPACK_dtrsyl("N", "N", &isgn, &m, &n, A.base(), &m, B.base(), &n,
lapack_int info;
dtrsyl("N", "N", &isgn, &m, &n, A.base(), &m, B.base(), &n,
C.base(), &m, &scale, &info);
if (info < -1)
throw SYLV_MES_EXCEPTION("Wrong parameter to LAPACK dtrsyl.");
......
......@@ -5,8 +5,8 @@
#include "SylvException.h"
#include "SylvMatrix.h"
#include "cppblas.h"
#include "cpplapack.h"
#include <dynblas.h>
#include <dynlapack.h>
#include <stdio.h>
#include <string.h>
......@@ -44,15 +44,15 @@ void SylvMatrix::multLeft(int zero_cols, const GeneralMatrix& a, const GeneralMa
// another copy of (usually big) b (we are not able to do inplace
// submatrix of const GeneralMatrix)
if (a.getLD() > 0 && ld > 0) {
int mm = a.numRows();
int nn = cols;
int kk = a.numCols();
blas_int mm = a.numRows();
blas_int nn = cols;
blas_int kk = a.numCols();
double alpha = 1.0;
int lda = a.getLD();
int ldb = ld;
blas_int lda = a.getLD();
blas_int ldb = ld;
double beta = 0.0;
int ldc = ld;
BLAS_dgemm("N", "N", &mm, &nn, &kk, &alpha, a.getData().base(), &lda,
blas_int ldc = ld;
dgemm("N", "N", &mm, &nn, &kk, &alpha, a.getData().base(), &lda,
b.getData().base()+off, &ldb, &beta, data.base(), &ldc);
}
}
......@@ -206,29 +206,30 @@ void SqSylvMatrix::multInvLeft2(GeneralMatrix& a, GeneralMatrix& b,
}
// PLU factorization
Vector inv(data);
int * const ipiv = new int[rows];
int info;
LAPACK_dgetrf(&rows, &rows, inv.base(), &rows, ipiv, &info);
lapack_int * const ipiv = new lapack_int[rows];
lapack_int info;
lapack_int rows2 = rows;
dgetrf(&rows2, &rows2, inv.base(), &rows2, ipiv, &info);
// solve a
int acols = a.numCols();
lapack_int acols = a.numCols();
double* abase = a.base();
LAPACK_dgetrs("N", &rows, &acols, inv.base(), &rows, ipiv,
abase, &rows, &info);
dgetrs("N", &rows2, &acols, inv.base(), &rows2, ipiv,
abase, &rows2, &info);
// solve b
int bcols = b.numCols();
lapack_int bcols = b.numCols();
double* bbase = b.base();
LAPACK_dgetrs("N", &rows, &bcols, inv.base(), &rows, ipiv,
bbase, &rows, &info);
dgetrs("N", &rows2, &bcols, inv.base(), &rows2, ipiv,
bbase, &rows2, &info);
delete [] ipiv;
// condition numbers
double* const work = new double[4*rows];
int* const iwork = new int[rows];
lapack_int* const iwork = new lapack_int[rows];
double norm1 = getNorm1();
LAPACK_dgecon("1", &rows, inv.base(), &rows, &norm1, &rcond1,
dgecon("1", &rows2, inv.base(), &rows2, &norm1, &rcond1,
work, iwork, &info);
double norminf = getNormInf();
LAPACK_dgecon("I", &rows, inv.base(), &rows, &norminf, &rcondinf,
dgecon("I", &rows2, inv.base(), &rows2, &norminf, &rcondinf,
work, iwork, &info);
delete [] iwork;
delete [] work;
......
......@@ -6,8 +6,8 @@
#include "SylvException.h"
#include "KronVector.h"
#ifdef MATLAB
#include "mex.h"
#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
# include <dynmex.h>
#endif
#include <math.h>
......@@ -30,7 +30,7 @@ void SylvMemoryPool::init(size_t size)
#ifdef USE_MEMORY_POOL
length = size;
#ifdef MATLAB
#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
if (base)
throw SYLV_MES_EXCEPTION("Attempt to use matlab memory pool twice.");
base = (char*) mxMalloc(length);
......@@ -90,7 +90,7 @@ SylvMemoryPool::~SylvMemoryPool()
void SylvMemoryPool::reset()
{
#ifndef MATLAB
#if !defined(MATLAB_MEX_FILE) && !defined(OCTAVE_MEX_FILE)
delete [] base;
base = 0;
allocated = 0;
......@@ -134,7 +134,7 @@ void operator delete[](void* p)
#ifdef USE_MEMORY_POOL
void* MallocAllocator::operator new(size_t size)
{
#ifdef MATLAB
#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
throw SYLV_MES_EXCEPTION("Attempt to call wrong memory allocator.");
#else
void* res = malloc(size);
......@@ -146,7 +146,7 @@ void* MallocAllocator::operator new(size_t size)
void* MallocAllocator::operator new[](size_t size)
{
#ifdef MATLAB
#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
throw SYLV_MES_EXCEPTION("Attempt to call wrong memory allocator.");
#else
void* res = malloc(size);
......@@ -158,7 +158,7 @@ void* MallocAllocator::operator new[](size_t size)
void MallocAllocator::operator delete(void* p)
{
#ifdef MATLAB
#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
throw SYLV_MES_EXCEPTION("Attempt to call wrong memory destructor.");
#else
free(p);
......@@ -167,7 +167,7 @@ void MallocAllocator::operator delete(void* p)
void MallocAllocator::operator delete[](void* p)
{
#ifdef MATLAB
#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
throw SYLV_MES_EXCEPTION("Attempt to call wrong memory destructor.");
#else
free(p);
......
......@@ -130,7 +130,7 @@ void SylvParams::setArrayNames(int& num, const char** names) const
names[num++] = "cpu_time";
}
#ifdef MATLAB
#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
mxArray* SylvParams::DoubleParamItem::createMatlabArray() const
{
return mxCreateScalarDouble(value);
......@@ -164,7 +164,7 @@ mxArray* SylvParams::createStructArray() const
const char* names[50];
int num;
setArrayNames(num, names);
const int dims[] = {1, 1};
const mwSize dims[] = {1, 1};
mxArray* const res = mxCreateStructArray(2, dims, num, names);
int i = 0;
......
......@@ -8,8 +8,8 @@
#include <stdio.h>
#include <string.h>
#ifdef MATLAB
#include "mex.h"
#if defined(MATLAB_MEX_FILE) || defined(OCTAVE_MEX_FILE)
# include <dynmex.h>
#endif
typedef enum {def, changed, undef} status;
......@@ -63,7 +63,7 @@ protected:
DoubleParamItem(const DoubleParamItem& item) : ParamItem<double>(item) {}
const DoubleParamItem& operator=(const double& val