Skip to content

Instantly share code, notes, and snippets.

/dmatrixv.cxx

Created February 8, 2014 01:09
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save anonymous/8875125 to your computer and use it in GitHub Desktop.
Save anonymous/8875125 to your computer and use it in GitHub Desktop.
Download ZIP
Raw
dmatrixv.cxx
/*********************************************************************************************
This file is part of the DMatrix library, a C++ tool for numerical linear algebra
Copyright (C) 2009 Victor M. Becerra
This library is free software; you can redistribute it and/or
modify it under the terms of the GNU Lesser General Public
License as published by the Free Software Foundation; either
version 2.1 of the License, or (at your option) any later version.
This library 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
Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public
License along with this library; if not, write to the Free Software
Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA,
or visit http://www.gnu.org/licenses/
Author: Dr. Victor M. Becerra
University of Reading
School of Systems Engineering
P.O. Box 225, Reading RG6 6AY
United Kingdom
e-mail: v.m.becerra@reading.ac.uk
**********************************************************************************************/
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <time.h>
#include <cstdio>
#include <stdlib.h>
#include "dmatrixv.h"
#define ERROR_MESSAGE error_message
#define PRINTF printf
#ifndef MAX
#define MAX(a, b) ( (a)>(b)? (a):(b) )
#endif
#ifndef MIN
#define MIN(a, b) ( (a)<(b)? (a):(b) )
#endif
DMatrix* DMatrix::GetTempPr(int i) {
auxPr[i].mtype = 0;
auxPr[i].mt = NULL;
auxPr[i].rowIndx = NULL;
auxPr[i].colIndx = NULL;
return &auxPr[i];
}
inline long ChkTmpIndx( long taindx )
{
if ( taindx >= N_TEMP_OBJECTS )
ERROR_MESSAGE(" Temporary arrays error: Increase N_TEMP_OBJECTS");
#ifdef DEBUG_TEMPS
printf("\n Temp Created --> indx: %d", taindx );
#endif
return taindx;
}
/* Definition of static member variables of class DMatrix */
DMatrix* DMatrix::auxPr = NULL; // Pointer to auxiliary arrays
int DMatrix::initFlag = 0; // To be allocated in AllocateAuxArr() */
int DMatrix::noAuxArr=N_TEMP_OBJECTS;// Number of auxiliary arrays
long DMatrix::dimAux=D_TEMP_OBJECTS; // Dimension of auxiliary arrays
int DMatrix::auxIndx=-1; // Index of used auxiliary arrays
int DMatrix::memberFlag = 0; // Member function flag
int DMatrix::errorFlag = 0; // Error flag
int DMatrix::print_level = 1; // Print level
const double DMatrix::MACH_EPS = MachEps();
time_t DMatrix::start_time = time(NULL);
clock_t DMatrix::start_clock = 0;
// long DMatrix::seed[RAND_STREAMS] = {RAND_DEFAULT};
long generate_seed(time_t *timer);
long DMatrix::seed[RAND_STREAMS] = {generate_seed(&DMatrix::start_time)};
int DMatrix::stream = 0; // stream index
#ifdef OTISS_CASE
double DMatrix::axp[N_TEMP_OBJECTS][D_TEMP_OBJECTS];
#endif
char* num2str(double num);
long generate_seed(time_t *timer)
{
long seed;
tm* st = localtime(timer);
seed = st->tm_yday*365*24*3600;
seed += st->tm_hour*3600;
seed += st->tm_min*60;
seed += st->tm_sec;
seed *=37;
return seed;
}
/* Declaration of Dummy object to allocate and de-allocate */
/* auxiliary arrays */
#ifndef MATLAB_MEX_FILE
#ifndef DECLARED_TEMPS
InitializeDMatrixClass Dummy;
#endif
#endif
#define NR_END 1
#define FREE_ARG char*
DMatrix::DMatrix(void)
// Default constructor
{
n = 0;
m = 0;
a = NULL;
asize = 0;
atype = 0;
auxFlag=0;
SetReferencedDMatrixPointer( NULL );
SetRowIndexPointer(NULL);
SetColIndexPointer(NULL);
SetMType(0);
allocated = false;
}
DMatrix::DMatrix(long Initn, long Initm)
// Matrix constructor using memory allocation
{
n=Initn;
m=Initm;
asize = n*m;
if ( Initm<0 || Initn<0 ) {
ERROR_MESSAGE("Attempt to create a matrix with \
negative dimension in \
DMatrix::DMatrix(long,long)");
}
if (Initm!=0 && Initn!=0) {
DMatrix::Allocate(asize);
memset(a, 0, asize*sizeof(double) );
}
atype = 0;
auxFlag = 0;
SetReferencedDMatrixPointer( NULL );
SetRowIndexPointer(NULL);
SetColIndexPointer(NULL);
SetMType(0);
}
DMatrix::DMatrix( const DMatrix& A)
// copy constructor
{
n=A.n;
m=A.m;
asize = n*m;
DMatrix::Allocate(asize);
memcpy( a, A.a , n*m*sizeof(double) );
atype = A.atype;
SetReferencedDMatrixPointer( NULL );
SetRowIndexPointer(NULL);
SetColIndexPointer(NULL);
SetMType(0);
auxFlag=0;
if( !DMatrix::GetMemberFlag() ) DMatrix::SetAuxIndx( -1 );
}
void DMatrix::PrintInfo(const char *str) const
{
fprintf(stderr,"\nInformation about DMatrix object: %s", str);
fprintf(stderr,"\nNumber of Rows :\t\t%ld",GetNoRows() );
fprintf(stderr,"\nNumber of Columns :\t\t%ld",GetNoCols() );
fprintf(stderr,"\nAllocation type :\t\t%d",atype);
fprintf(stderr,"\nAllocation size :\t\t%ld",asize);
}
void DMatrix::AllocateAuxArr( void )
{
// Routine to allocate the temporary objects at program start
int i;
if (DMatrix::GetInitFlag()==0 ) {
*(DMatrix::GetAuxPr()) = (DMatrix*) my_calloc( N_TEMP_OBJECTS , sizeof(DMatrix) );
for (i =0; i< GetNoAuxArr(); i++ )
{
#ifdef DECLARED_TEMPS
DMatrix::auxPr[i].a = axp[i];
DMatrix::auxPr[i].atype = 1;
#endif /* DECLARED_TEMPS */
DMatrix::auxPr[i].Resize(GetDimAux(), 1 );
DMatrix::auxPr[i].auxFlag = 1;
DMatrix::auxPr[i].asize = GetDimAux();
}
DMatrix::SetInitFlag( 1 );
} // End if
}
void DMatrix::DeAllocateAuxArr( void )
{
// Routine to deallocate temporary objects at program termination
int i;
for ( i=0; i< GetNoAuxArr(); i++ )
{
DMatrix::auxPr[i].~DMatrix();
}
mxFree( *(DMatrix::GetAuxPr()) );
}
void DMatrix::Resize( long nnrow, long nncol )
{
// Resizes the row and column dimensions of a DMatrix object
// additional memory is allocated if necessary
long i;
if (nnrow==0 || nncol==0) {
n = nnrow;
m = nncol;
if (atype==0 && auxFlag!=1) {
DMatrix::DeAllocate();
asize = 0;
}
return;
}
if (n==nnrow && m==nncol ) return;
if ( n*m == nnrow*nncol)
{
n = nnrow;
m = nncol;
return;
}
if (nnrow < 0 || nncol < 0 ) { error_message("Negative arguments in Resize()"); }
if (asize >= nnrow*nncol && atype==1 )
{
for (i=nnrow*nncol+1;i<=asize;i++) {
a[i-1]=0.0;
}
n=nnrow;
m=nncol;
return;
}
if ( a!=NULL && atype!=1 && auxFlag!=1 )
{
double* atemp = (double *) my_calloc( n*m, sizeof(double) );
if (!atemp) {
ERROR_MESSAGE("\nError allocating double array in DMatrix::Resize()");
}
memcpy( atemp, a, n*m*sizeof(double) );
DMatrix::DeAllocate();
asize = nnrow*nncol;
DMatrix::Allocate(asize);
memcpy( a, atemp, MIN( n*m, asize )*sizeof(double) );
n = nnrow;
m = nncol;
mxFree( atemp ) ;
return;
}
if ( a!=NULL && auxFlag==1 ) {
if (nnrow*nncol <= dimAux )
{
n = nnrow;
m = nncol;
}
else {
ERROR_MESSAGE("\nError resizing Temporary array in DMatrix::Resize()");
}
}
if ( a==NULL && atype==0 )
{
n = nnrow;
m = nncol;
asize = nnrow*nncol;
DMatrix::Allocate(asize);
return;
}
if ( atype == 1 && nnrow*nncol > asize )
{
ERROR_MESSAGE("Size error in DMatrix::Resize");
}
}
void DMatrix::SetPrintLevel( int plevel )
{
print_level=plevel;
return;
}
int DMatrix::PrintLevel() {
return print_level;
}
void DMatrix::Print(const char* text) const
{
// Prints a matrix
// Eg. A.Print(" A = ");
long i,j;
if (print_level) {
#ifndef MATLAB_MEX_FILE
fprintf(OUTPUT_STREAM,"\n%s\n",text);
for (i=1;i<=n;i++)
{
for (j=1;j<=m;j++)
{
fprintf(OUTPUT_STREAM," %s ", num2str(elem(i,j)) );
fprintf(OUTPUT_STREAM,"\t");
}
fprintf(OUTPUT_STREAM,"\n");
}
#else
mexPrintf("\n%s\n",text);
for (i=1;i<=n;i++)
{
for (j=1;j<=m;j++)
{
mexPrintf(" %s ", num2str(elem(i,j)) );
mexPrintf("\t");
}
mexPrintf("\n");
}
#endif
}
if( !DMatrix::GetMemberFlag() ) DMatrix::SetAuxIndx( -1 );
}
void error_message(const char *error_text)
{
// Error handling routine
string m1;
string m2;
m1 = "\n**** Run time error";
m1 += "\n**** To trace this error, set up your debugger to break at";
m1 += "\n**** function 'error_message', then do a backtrace.";
m1 += "\n**** A diagnostic message is given below:";
m2 = error_text;
m2 = "\n**** ====> " + m2 + " <====\n\n";
if ( DMatrix::PrintLevel() ) {
#ifndef MATLAB_MEX_FILE
fprintf(OUTPUT_STREAM,"%s", m1.c_str());
fprintf(OUTPUT_STREAM,"%s", m2.c_str());
FILE* err_file = fopen("error_message.txt","w");
fprintf(err_file,"%s", m1.c_str() );
fprintf(err_file,"%s", m2.c_str() );
fclose(err_file);
#else
mexPrintf("\n%s\n",m1.c_str());
mexPrintf("%s\n",m2.c_str());
#endif
}
DMatrix::RiseErrorFlag();
throw ErrorHandler(m1+m2);
}
ErrorHandler::ErrorHandler(const string m)
{
error_message = m;
}
DMatrix& DMatrix::operator- (const DMatrix& Other_matrix) const
{
// Matrix substraction operator: substracts two matrices
// Example C = D - E; where C,D,E are DMatrix objects
long i;
long nelem = n*m;
const double *OtherPr = Other_matrix.GetConstPr();
DMatrix* Temp;
if (m != Other_matrix.m || n != Other_matrix.n)
{ ERROR_MESSAGE("Incoherent matrix dimensions in substraction"); }
Temp = DMatrix::GetTempPr(ChkTmpIndx(DMatrix::IncrementAuxIndx()));
Temp->Resize( n, m );
for (i=0; i< nelem; i++)
{
Temp->GetPr()[i] = a[i] - OtherPr[i];
}
return *Temp;
}
DMatrix& DMatrix::operator/ (double Arg) const
{
// Matrix scalar division.
// Eg. C = A/x, where C & A are DMatrix objects and x is a double
long i;
long nelem = n*m;
DMatrix* Temp;
double *pr;
Temp = DMatrix::GetTempPr(ChkTmpIndx(DMatrix::IncrementAuxIndx()));
Temp->Resize( n, m );
pr = Temp->GetPr();
for( i=0; i< nelem; i++ )
{
pr[i] = a[i]/Arg;
}
return *Temp;
}
DMatrix& operator *(double Arg, const DMatrix& A)
{
// Scalar by Matrix product.
// Eg. C = x*A, where C & A are DMatrix objects and x is a double
long i,n,m,nelem;
const double * Apr = A.GetConstPr();
n=A.GetNoRows();
m=A.GetNoCols();
nelem = n*m;
DMatrix* Temp;
double *pr;
Temp = DMatrix::GetTempPr(ChkTmpIndx(DMatrix::IncrementAuxIndx()));
Temp->Resize( n, m );
pr = Temp->GetPr();
for (i=0; i< nelem; i++ )
{
pr[i] = Apr[i]*Arg;
}
return *Temp;
}
void DMatrix::SetMType( int arg )
{
if (arg !=0 && arg !=1) {
ERROR_MESSAGE("Incorrect argument value in DMatrix::SetMType");
}
mtype = arg;
}
DMatrix& DMatrix::operator= (const DMatrix& Other_matrix)
{
// Matrix Assignment operator
// Eg. A = B;
if ( GetMType() == 1 ) {
return (AssignmentToColonReference( Other_matrix ));
}
else {
if (n > 0 && m>0) {
if (n != Other_matrix.n || m != Other_matrix.m)
this->Resize(Other_matrix.n,Other_matrix.m);
}
else {
n = Other_matrix.n;
m = Other_matrix.m;
DMatrix::Allocate(n*m);
}
memcpy( a, Other_matrix.a, n*m*sizeof(double) );
if( !DMatrix::GetMemberFlag() ) DMatrix::SetAuxIndx( -1 );
return *this;
} // End if //
}
void DMatrix::DeAllocate()
// free the arrays allocated with Allocate() //
{
if (n!=0 && m!=0 ) {
if (allocated==true) {
mxFree((FREE_ARG) (a) );
allocated=false;
}
}
}
DMatrix& zeros(long nrows, long ncols)
{
long nelem = nrows*ncols;
double *pr;
DMatrix* Temp;
Temp = DMatrix::GetTempPr(ChkTmpIndx(DMatrix::IncrementAuxIndx()));
Temp->Resize( nrows, ncols );
pr = Temp->GetPr();
memset(pr, 0, nelem*sizeof(double) );
return *Temp;
}
DMatrix& ones(long nrows, long ncols)
{
long i;
long nelem = nrows*ncols;
double *pr;
DMatrix* Temp;
Temp = DMatrix::GetTempPr(ChkTmpIndx(DMatrix::IncrementAuxIndx()));
Temp->Resize( nrows, ncols );
pr = Temp->GetPr();
for ( i=0; i< nelem; i++ ) {
pr[i] = 1.0;
}
return *Temp;
}
DMatrix& DMatrix::operator &(const DMatrix B) const
{
return elemProduct(*this,B);
}
DMatrix& elemProduct( const DMatrix& A, const DMatrix& B )
{
// element by element matrix product
DMatrix* Temp;
Temp = DMatrix::GetTempPr(ChkTmpIndx(DMatrix::IncrementAuxIndx()));
Temp->Resize( A.n, A.m );
if ( (A.n != B.n) || (A.m != B.m) ) {
error_message( "Dimension error in elemProduct()");
}
for ( long i=0; i< A.n*A.m; i++ ) {
Temp->a[i] = A.a[i] * B.a[i] ;
}
return *Temp;
}
double& DMatrix::operator() (long row, long col )
{
if (row > n || col > m )
DMatrix::Resize(row,col);
if (row < 0 || col < 0)
ERROR_MESSAGE("negative index error in A(row,col)");
return elem( row, col );
}
double& DMatrix::operator() (long index )
{
if (index > asize && (n==1 || m==1) )
{
if (n==1)
{
DMatrix::Resize(1,index);
}
else if (m==1)
{
DMatrix::Resize(index,1);
}
}
if ( (n==0 || m==0) )
{
DMatrix::Resize(index,1);
}
if (index > asize && (n>1 && m>1) )
{
ERROR_MESSAGE("Non-vector matrix cannot be resized in call to A(I) in LHS");
}
if (index > asize || index < 1 ) ERROR_MESSAGE("Index error in operator()(long)");
return a[ index - 1 ];
}
double DMatrix::operator() (long index ) const
{
if (index > asize || index < 1 ) ERROR_MESSAGE("Index error in operator()(long)");
return a[ index - 1 ];
}
DMatrix::~DMatrix()
{
if (atype==0 && a!=NULL )
{
DMatrix::DeAllocate();
}
}
double Max(const DMatrix & A, int* rindx, int* cindx )
{
double mx;
const double *Apr = A.GetConstPr();
int ri, ci;
long inx;
mx = Apr[0];
inx = 0;
for (long i=1; i<A.n*A.m; i++)
{
if (Apr[i] > mx) { mx = Apr[i]; inx = i; }
}
ci = (inx)/A.n + 1;
ri = (inx+1) - (ci-1)*A.n;
if ( rindx != NULL ) *rindx = ri;
if ( cindx != NULL ) *cindx = ci;
return mx;
}
void DMatrix::Allocate(long size)
{
/* Memory allocation function */
a = (double *) my_calloc( size, sizeof(double) );
if (!a) ERROR_MESSAGE("allocation failure in DMatrix::Allocate()");
// a = a-nl+NR_END;
asize = size;
allocated = true;
}
// extern "C" char *gcvt( double, int, char*);
char* num2str(double num)
{
static char str[20];
#ifndef MEX_OR_XPC
int sig = 10;
gcvt( num, sig, str);
#endif
return str;
}
double MachEps(void)
{
// Function to compute machine precision
double meps = 1.,eps = 2.;
while(eps != 1.)
{
meps /= 2.;
eps = 1. + meps;
}
return 2.0*meps;
}
// Indexing functions
DMatrix& DMatrix::operator() (const DMatrix& RowIndx, const DMatrix& ColIndx )
{
long nrows;
long ncols;
long ii;
long jj;
int rflag = 0;
int cflag = 0;
DMatrix* RtVal;
if ( RowIndx(1) == -1 ) rflag = 1;
if ( ColIndx(1) == -1 ) cflag = 1;
if( !rflag ) nrows = RowIndx.GetNoRows();
else nrows = GetNoRows();
if( !cflag ) ncols = ColIndx.GetNoRows();
else ncols = GetNoCols();
RtVal = DMatrix::GetTempPr(ChkTmpIndx(DMatrix::IncrementAuxIndx()));
RtVal->SetMType( 1 );
RtVal->SetReferencedDMatrixPointer( this );
RtVal->SetRowIndexPointer( &RowIndx );
RtVal->SetColIndexPointer( &ColIndx );
RtVal->Resize( nrows, ncols );
for (long i = 1; i<=nrows; i++ ) {
if( !rflag ) ii = (long) RowIndx(i);
else ii = i;
for (long j = 1; j<=ncols; j++ ) {
if (!cflag) jj = (long) ColIndx(j);
else jj=j;
RtVal->elem(i,j) = elem(ii,jj);
}
}
return (*RtVal);
}
DMatrix& DMatrix::operator() (long row, const DMatrix& ColIndx )
{
DMatrix* RowIndx;
RowIndx = DMatrix::GetTempPr(ChkTmpIndx(DMatrix::IncrementAuxIndx()));
RowIndx->SetMType( 0 );
RowIndx->Resize( 1 , 1 );
RowIndx->elem( 1,1 ) = (double) row;
return (this->operator()(*RowIndx,ColIndx) );
}
DMatrix& DMatrix::AssignmentToColonReference( const DMatrix& RightMt )
{
long ii,jj;
int rflag = 0;
int cflag = 0;
if ( (*rowIndx)(1) == -1.0 ) rflag = 1;
if (colIndx!=NULL) {
if ( (*colIndx)(1) == -1.0 ) cflag = 1;
}
if (GetNoRows() != RightMt.GetNoRows() ||
GetNoCols() != RightMt.GetNoCols() )
{
ERROR_MESSAGE("Dimensions error in AssignmentToColonReference()");
}
if ( Max(*rowIndx)>mt->GetNoRows() )
{
mt->Resize( (long) Max(*rowIndx) , mt->GetNoCols() );
}
if (colIndx!=NULL) {
if ( Max(*colIndx)>mt->GetNoCols() )
{
mt->Resize( mt->GetNoRows() , (long) Max(*colIndx) );
}
}
for (long i = 1; i<=GetNoRows(); i++ ) {
if (!rflag) ii = (long) (*rowIndx)(i);
else ii = i;
for (long j = 1; j<=GetNoCols(); j++ ) {
if (colIndx!=NULL) {
if (!cflag) jj = (long) (*colIndx)(j);
else jj = j;
} else {jj=1;}
mt->elem(ii,jj) = RightMt.elem(i,j);
}
}
if( !DMatrix::GetMemberFlag() ) DMatrix::SetAuxIndx( -1 );
return (*this);
}
DMatrix& colon( void )
{
DMatrix* Temp;
Temp = DMatrix::GetTempPr(ChkTmpIndx(DMatrix::IncrementAuxIndx()));
Temp->Resize( 1, 1 );
Temp->elem(1,1) = -1.0;
return (*Temp);
}
long length(const DMatrix& A)
{
int retval;
retval = A.GetNoRows()*A.GetNoCols();
return retval;
}
void* my_calloc(size_t num, size_t size )
{
void* pointer;
pointer = mxCalloc(num, size);
if (pointer==NULL) {
error_message("Memory allocation error in my_calloc()");
}
return pointer;
}
Raw
dmatrixv.h
/************************************************************************************************
This file is part of the DMatrix library, a C++ tool for numerical linear algebra
Copyright (C) 2009 Victor M. Becerra
This library is free software; you can redistribute it and/or
modify it under the terms of the GNU Lesser General Public
License as published by the Free Software Foundation; either
version 2.1 of the License, or (at your option) any later version.
This library 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
Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public
License along with this library; if not, write to the Free Software
Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA,
or visit http://www.gnu.org/licenses/
Author: Dr. Victor M. Becerra
University of Reading
School of Systems Engineering
P.O. Box 225, Reading RG6 6AY
United Kingdom
e-mail: v.m.becerra@ieee.org
**********************************************************************************************/
/*! \mainpage DMatrix and SparseMatrix classes
\section intro Introduction
The author developed the main features of the DMatrix class between 1994 and 1999. The class makes extensive use of operator overloading in order to facilitate the implementation in C++ of complicated matrix expressions, and it has intefaces to a number of LAPACK routines. In 2008, the class has been tested with current compilers, its functionality was expanded, and the code was published under the GNU Lesser General Public License. The class at present is restricted to dense and real matrices. In 2008, the SparseMatrix class was added to the library to incorporate basic sparse matrix functionality. The SparseMatrix class offers interfaces to some functions available in the CXSparse and LUSOL libraries.
The library should compile without problems with the following C++ compilers: GNU C++ version 4.X and Microsoft Visual Studio 2005.
\section license License
This work is copyright (c) Victor M. Becerra (2009)
This library is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public
License as published by the Free Software Foundation; either version 2.1 of the License, or (at your option) any later version.
This library 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 Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public License along with this library; if not, write to the Free Software
Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA, or visit http://www.gnu.org/licenses/
Author: Dr. Victor M. Becerra, University of Reading, School of Systems Engineering, P.O. Box 225, Reading RG6 6AY, United Kingdom, e-mail: v.m.becerra@ieee.org.
\section install Installing the library
See the INSTALL file.
\section examples Examples of use
See the source code in the examples directory.
*/
#ifndef DMATRIX_H
#define DMATRIX_H
#ifndef bool
// typedef int bool;
#endif
#ifndef N_TEMP_OBJECTS
#define N_TEMP_OBJECTS (40)
#endif
#ifndef D_TEMP_OBJECTS
#define D_TEMP_OBJECTS (1000000)
#endif
#ifndef OUTPUT_STREAM
#define OUTPUT_STREAM stderr
#endif
#include <stdio.h>
#include <stdarg.h>
#ifndef MATLAB_MEX_FILE
#define mxCalloc calloc
#define mxFree free
#endif /* MATLAB_MEX_FILE */
#include <time.h>
#ifndef true
#define true 1
#endif
#ifndef false
#define false 0
#endif
#include <string>
using std::string;
#define RAND_MODULUS 2147483647 /* DO NOT CHANGE THIS VALUE */
#define RAND_MULTIPLIER 48271 /* DO NOT CHANGE THIS VALUE */
#define RAND_CHECK 399268537 /* DO NOT CHANGE THIS VALUE */
#define RAND_STREAMS 256 /* # of streams, DO NOT CHANGE THIS VALUE */
#define RAND_A256 22925 /* jump multiplier, DO NOT CHANGE THIS VALUE */
#define RAND_DEFAULT 123456789 /* initial seed, use 0 < RAND_DEFAULT < RAND_MODULUS */
class DMatrix;
class SparseMatrix;
//! DMatrix class
/**
A C++ class for dense and real matrix and vector computations with interfaces to
a number of LAPACK functions
*/
class DMatrix { // define DMatrix class
protected:
//! Array of doubles to store matrix elements using column major storage
double *a;
//! Number of matrix rows
long n;
//! Number of matrix columns
long m;
//! Number of allocated elements in a
long asize;
/*! Flag to indicate type of allocation.
type = 0 : allocated matrix
type = 1 : non-allocated matrix, uses predefined array for storage
*/
int atype;
/*! Flag to indicate type of matrix
mtype = 0 : normal matrix
mtype = 1 : colon - reference matrix
*/
int mtype;
//! Flag to indicate auxiliary (temporary) matrix flag = 1 if temporary matrix, 0 otherwise.
int auxFlag;
//! Flag to indicate that element storage has been allocated
int allocated;
//! Referenced matrix pointer
DMatrix* mt;
//! Row indices
const DMatrix* rowIndx;
//! Column indices
const DMatrix* colIndx;
// Protected static members
//! Array of Temporary matrices
static DMatrix* auxPr;
//! Number of auxiliary matrices
static int noAuxArr;
//! Dimension of auxiliary arrays
static long dimAux;
//! Index of used auxiliary matrices
static int auxIndx;
//! Member function flag to control resetting of auxIndx
static int memberFlag;
//! Flag to indicate aux arrays allocation
static int initFlag;
//! Machine precision constant
static const double MACH_EPS;
//! Flag to indicate error condition
static int errorFlag;
//! Print level flag, 1: output sent to sderr, 0: no output sent
static int print_level;
//! current state of each stream
static long seed[];
//! stream index for pseudo-ramdon number generator
static int stream;
//! variable to store start time after tic() call.
static time_t start_time;
//! clock_t variable
static clock_t start_clock;
//! returns pointer to the i-th temporary object
static DMatrix* GetTempPr(int i);
//! Gets the memberFlag value from the object
static int GetMemberFlag() { return memberFlag; }
//! Sets the memberFlag value
static void SetMemberFlag(int arg) { memberFlag=arg; }
//! Checks the auxiliary objects
static void ChkAuxArrays();
//! Gets the number of temporary objects
static int GetNoAuxArr() { return noAuxArr; }
//! Gets the dimensions of the array of temporary objects
static long GetDimAux() { return dimAux; }
//! Gets the value of initFlag
static int GetInitFlag() { return initFlag; }
//! Sets the value of initFlag
static void SetInitFlag(int arg) { initFlag = arg ; }
//! Gets the current index of temporary objects
static int GetAuxIndx() { return auxIndx; }
//! Increments the index of temporary objects
static int IncrementAuxIndx() {
#ifdef CHECK_AUX
fprintf(stderr,"\nauxIndx=%d",auxIndx);
#endif
return ++auxIndx; }
//! Decrements the index of temporary objects
static int DecrementAuxIndx() { return --auxIndx ; }
//! Sets the index of temporary objects
static void SetAuxIndx( int i ) { auxIndx = i; }
//! Sets the dimension of each object in the array of temporary objects
static void SetDimAux( long dd ) { dimAux = dd; }
//! Sets the number of auxiliary objects
static void SetNoAuxArr( int nn ) { noAuxArr = nn ; }
//! Sets the errorFlag member to true
static void RiseErrorFlag() { errorFlag = true; }
//! Sets the value of auxFlag
void SetAuxFlag(int arg) { auxFlag = arg; }
//! Gets the value of auxFlag
int GetAuxFlag() { return auxFlag; }
//! Gets the value of start_clock member
static clock_t GetStartTicks(void) { return start_clock; }
//! Sets the value of start_clock member
static void SetStartTicks(clock_t st) { start_clock=st; }
#ifdef DECLARED_TEMPS
static double axp[N_TEMP_OBJECTS][D_TEMP_OBJECTS];
#endif
#define MC_EPSILON 2.221e-16
void Allocate(long size);
void DeAllocate(); // De-allocate memory
void SetReferencedDMatrixPointer( DMatrix* arg ) { mt = arg; }
DMatrix* GetReferencedDMatrixPointer() { return mt; }
void SetRowIndexPointer(const DMatrix* arg ) { rowIndx = arg; }
void SetColIndexPointer(const DMatrix* arg ) { colIndx = arg; }
void SetMType( int arg );
int GetMType() { return mtype; }
DMatrix& AssignmentToColonReference( const DMatrix& A );
DMatrix& AssignmentToColonReference( double arg );
public:
static void AllocateAuxArr( void );
static void DeAllocateAuxArr( void );
static double random_uniform(void);
static DMatrix** GetAuxPr(void) { return &auxPr; }
static int isThereError(void) { return errorFlag; }
void Print(const char *text) const;
void PrintInfo(const char *text) const;
static void SetPrintLevel( int plevel );
static int PrintLevel();
void FillWithZeros( void );
long getn() {return n;}
long getm() {return m;}
long getn() const {return n;}
long getm() const {return m;}
long GetNoRows() { return n; }
long GetNoRows() const { return n; }
long GetNoCols() { return m; }
long GetNoCols() const { return m; }
double *GetPr() { return a; }
double *GetConstPr() const { return a; }
double* geta() {return a;}
int getatype() { return atype; }
int isEmpty() { if (n*m == 0) return 1; else return 0; }
int isVector() const { if (n==1 || m==1) return 1; else return 0; }
double element(long i, long j );
double& elem( long i, long j ) { return a[ (j-1)*n + i-1 ]; }
double elem( long i, long j ) const { return a[(j-1)*n+i-1]; }
DMatrix(void); // Default constructor
DMatrix( long Initn, long Initm);
DMatrix( long vDim, double* v, long Initn, long Initm );
DMatrix( long Initn); // Column vector constructor using memory allocation
DMatrix(long Initn,long Initm,double a11,...);
DMatrix( const DMatrix& A); // copy constructor
~DMatrix();
void Resize(long nnrow, long nncol );
DMatrix& operator+ (const DMatrix& rval) const;
DMatrix& operator += (const DMatrix &rval);
DMatrix& operator+ (double x) const;
DMatrix& operator- (const DMatrix& rval) const;
DMatrix& operator- (double x) const;
friend DMatrix& operator- (const DMatrix& A);
DMatrix& operator* (const DMatrix& rval) const;
DMatrix& operator* (double Arg) const;
DMatrix& operator/ (double Arg) const;
DMatrix& operator/ (const DMatrix& rval) const;
DMatrix& operator= (const DMatrix& rval);
DMatrix& operator = (double val);
DMatrix& operator = (const char* str);
DMatrix& operator &(const DMatrix B) const;
double& operator() (long row,long col);
double& operator() (long row,const char* end);
double& operator() (const char* end,long col);
double operator() (long row,long col) const;
double operator() (long row,const char* end) const;
double operator() (const char* end,long col) const;
double& operator() (long index );
double& operator() (const char* end);
double operator() (long k ) const;
double operator() (const char* end) const;
DMatrix& operator() (const DMatrix& RowIndx, const DMatrix& ColIndx );
DMatrix& operator() (const DMatrix& RowIndx );
DMatrix& operator() (const DMatrix& RowIndx, long col );
DMatrix& operator() (long row, const DMatrix& ColIndx);
DMatrix& operator() (const DMatrix& RowIndx,const char* end );
DMatrix& operator() (const char* end, const DMatrix& ColIndx);
friend DMatrix& colon( double i1, double increment, double i2 );
friend DMatrix& colon( int i1, int increment, int i2);
friend DMatrix& colon( int i1, int i2 );
friend DMatrix& colon( double i1, double i2);
friend DMatrix& colon( void );
friend DMatrix& operator *(double r, const DMatrix& A);
friend DMatrix& zeros(long n, long m);
friend DMatrix& ones(long n, long m);
friend double Max(const DMatrix& A,int* rindx=NULL, int* cindx=NULL);
friend DMatrix& elemProduct( const DMatrix& A, const DMatrix& B );
friend DMatrix& elemDivision( const DMatrix& A, const DMatrix& B);
static double GetEPS() { return MC_EPSILON ; }
friend void error_message(const char *input_text);
friend DMatrix& reshape(DMatrix& A, long N, long M);
friend long length(const DMatrix& A);
};
//! InitializeDMatrixClass class
/**
This is a dummy C++ class intended to initialise the temporary objects of the DMatrix class.
*/
class InitializeDMatrixClass{
public:
//! This is the default constructor which calls the function DMatrix::AllocateAuxArr().
InitializeDMatrixClass() { DMatrix::AllocateAuxArr(); }
//! This is the destructor which calls the function DMatrix::DeAllocateAuxArr().
~InitializeDMatrixClass() { DMatrix::DeAllocateAuxArr(); }
};
/* inline utility functions and prototypes */
// Swap double
inline void Swap(double *x,double *y)
{double temp = *x; *x = *y; *y = temp;}
// Swap int
inline void Swap(int *x,int *y)
{int temp = *x; *x = *y; *y = temp;}
double MachEps(void);
DMatrix& operator- (const DMatrix& A);
DMatrix& colon( double i1, double increment, double i2 );
DMatrix& colon( int i1, int increment, int i2);
DMatrix& colon( int i1, int i2 );
DMatrix& colon( double i1, double i2);
DMatrix& colon( void );
DMatrix& operator *(double Arg, const DMatrix& A);
DMatrix& zeros(long n, long m);
DMatrix& ones(long n, long m);
double Max(const DMatrix& A,int* rindx, int* cindx);
DMatrix& elemProduct( const DMatrix& A, const DMatrix& B );
DMatrix& elemDivision( const DMatrix& A, const DMatrix& B);
void error_message(const char *input_text);
DMatrix& reshape(DMatrix& A, long nn, long mm);
long length(const DMatrix& A);
void* my_calloc(size_t num, size_t size );
//! ErrorHandler class
/**
This is a C++ class intended to handle error conditions.
*/
class ErrorHandler
{
public:
//! A string of characters which contains the error message
string error_message;
//! A constructor which takes the error message as an argument and assigns it to error_message
/**
\param m is the error message string.
\sa function error_message().
*/
ErrorHandler(const string m);
};
#endif /* DMATRIX_H */
Raw
dmtest.cxx
#include "dmatrixv.h"
int main()
{
DMatrix x = zeros(1,20);
DMatrix pointx = zeros(1,20);
int i,j;
long n = length(pointx);
DMatrix L(n, x.GetNoCols() );
for(i=1;i<=n;i++) {
for(j=1;j<=x.GetNoCols();j++) {
L(i,j) = 1.0;
}
}
x(1, 1) = 0;
x(1, 2) = -9.628147553e-02;
x(1, 3) = -0.3203275059;
x(1, 4) = -0.6656101096;
x(1, 5) = -1.123158695;
x(1, 6) = -1.681117989;
x(1, 7) = -2.32503568;
x(1, 8) = -3.038234081;
x(1, 9) = -3.80224147;
x(1,10) = -4.597270314;
x(1,11) = -5.402729686;
x(1,12) = -6.19775853;
x(1,13) = -6.961765919;
x(1,14) = -7.67496432;
x(1,15) = -8.318882011;
x(1,16) = -8.876841305;
x(1,17) = -9.33438989;
x(1,18) = -9.679672494;
x(1,19) = -9.903718524;
x(1,20) = -10;
pointx(1, 1) = 0;
pointx(1, 2) = -0.5263157895;
pointx(1, 3) = -1.052631579;
pointx(1, 4) = -1.578947368;
pointx(1, 5) = -2.105263158;
pointx(1, 6) = -2.631578947;
pointx(1, 7) = -3.157894737;
pointx(1, 8) = -3.684210526;
pointx(1, 9) = -4.210526316;
pointx(1,10) = -4.736842105;
pointx(1,11) = -5.263157895;
pointx(1,12) = -5.789473684;
pointx(1,13) = -6.315789474;
pointx(1,14) = -6.842105263;
pointx(1,15) = -7.368421053;
pointx(1,16) = -7.894736842;
pointx(1,17) = -8.421052632;
pointx(1,18) = -8.947368421;
pointx(1,19) = -9.473684211;
pointx(1,20) = -10;
// code snippet from psopt.cxx function lagrange_interpolation
for (i=1;i<=n;i++) {
for (j=1;j<=n;j++) {
if (i != j) {
L(i,colon()) = ( L(i,colon())&( x-pointx(j)*ones(1,length(x)) ) )/(pointx(i)-pointx(j));
}
}
}
L.Print("");
return 0;
}
Raw
testscript.sh
#!/bin/sh
g++ -g -o dmtest-g++ -DUNIX dmatrixv.cxx dmtest.cxx
./dmtest-g++
clang++ -g -o dmtest-clang -DUNIX dmatrixv.cxx dmtest.cxx
./dmtest-clang
#lldb dmtest-clang
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment