/*! \file ccolmatrix.cpp
* \brief Compressed column sparse matrix algebra
*/
/* Copyright (c) 2005-2012 Taneli Kalvas. All rights reserved.
*
* You can redistribute this software and/or modify it under the terms
* of the GNU General Public License as published by the Free Software
* Foundation; either version 2 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
* General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this library (file "COPYING" included in the package);
* if not, write to the Free Software Foundation, Inc., 51 Franklin
* Street, Fifth Floor, Boston, MA 02110-1301 USA
*
* If you have questions about your rights to use or distribute this
* software, please contact Berkeley Lab's Technology Transfer
* Department at TTD@lbl.gov. Other questions, comments and bug
* reports should be sent directly to the author via email at
* taneli.kalvas@jyu.fi.
*
* NOTICE. This software was developed under partial funding from the
* U.S. Department of Energy. As such, the U.S. Government has been
* granted for itself and others acting on its behalf a paid-up,
* nonexclusive, irrevocable, worldwide license in the Software to
* reproduce, prepare derivative works, and perform publicly and
* display publicly. Beginning five (5) years after the date
* permission to assert copyright is obtained from the U.S. Department
* of Energy, and subject to any subsequent five (5) year renewals,
* the U.S. Government is granted for itself and others acting on its
* behalf a paid-up, nonexclusive, irrevocable, worldwide license in
* the Software to reproduce, prepare derivative works, distribute
* copies to the public, perform publicly and display publicly, and to
* permit others to do so.
*/
#include "config.hpp"
#include "ccolmatrix.hpp"
#include "crowmatrix.hpp"
#include "coordmatrix.hpp"
#include "mvector.hpp"
#include "sort.hpp"
#include <cstdlib>
#include <cmath>
#include <limits>
#include <iomanip>
inline void CColMatrix::allocate( void )
{
if( !(_ptr = (int *)malloc( (_m+1)*sizeof(int) )) ) {
_ptr = _row = NULL;
_val = NULL;
_n = _m = _nz = _asize = 0;
throw( ErrorNoMem( ERROR_LOCATION ) );
}
if( _asize == 0 ) {
_row = NULL;
_val = NULL;
} else {
if( !(_row = (int *)malloc( _asize*sizeof(int) )) ) {
free( _ptr );
_ptr = _row = NULL;
_val = NULL;
_n = _m = _nz = _asize = 0;
throw( ErrorNoMem( ERROR_LOCATION) );
}
if( !(_val = (double *)malloc( _asize*sizeof(double) )) ) {
free( _ptr );
free( _row );
_ptr = _row = NULL;
_val = NULL;
_n = _m = _nz = _asize = 0;
throw( ErrorNoMem( ERROR_LOCATION) );
}
}
}
inline void CColMatrix::reallocate( void )
{
int *tmp;
double *tmp2;
if( !(tmp = (int *)realloc( _ptr, (_m+1)*sizeof(int) )) ) {
free( _ptr );
free( _row );
free( _val );
_ptr = _row = NULL;
_val = NULL;
_n = _m = _nz = _asize = 0;
throw( ErrorNoMem( ERROR_LOCATION) );
}
_ptr = tmp;
if( _asize == 0 ) {
_row = NULL;
_val = NULL;
} else {
if( !(tmp = (int *)realloc( _row, _asize*sizeof(int) )) ) {
free( _ptr );
free( _row );
free( _val );
_ptr = _row = NULL;
_val = NULL;
_n = _m = _nz = _asize = 0;
throw( ErrorNoMem( ERROR_LOCATION) );
}
_row = tmp;
if( !(tmp2 = (double *)realloc( _val, _asize*sizeof(double) )) ) {
free( _ptr );
free( _row );
free( _val );
_ptr = _row = NULL;
_val = NULL;
_n = _m = _nz = _asize = 0;
throw( ErrorNoMem( ERROR_LOCATION) );
}
_val = tmp2;
}
}
CColMatrix::CColMatrix()
: _n(0), _m(0), _nz(0), _asize(0), _row(NULL), _val(NULL)
{
allocate();
}
CColMatrix::CColMatrix( int n, int m )
{
_n = n;
_m = m;
_nz = 0;
_asize = 0;
allocate();
memset( _ptr, 0, (_m+1)*sizeof(int) );
}
CColMatrix::CColMatrix( int n, int m, int nz, int asize,
int *ptr, int *row, double *val )
{
_n = n;
_m = m;
_nz = nz;
_asize = asize;
_ptr = ptr;
_row = row;
_val = val;
}
void CColMatrix::build( const class CColMatrix &mat )
{
_n = mat._n;
_m = mat._m;
_nz = mat._nz;
_asize = mat._nz;
reallocate();
memcpy( _ptr, mat._ptr, (_m+1)*sizeof(int) );
memcpy( _row, mat._row, _nz*sizeof(int) );
memcpy( _val, mat._val, _nz*sizeof(double) );
}
CColMatrix::CColMatrix( const CColMatrix &mat )
: _ptr(NULL), _row(NULL), _val(NULL)
{
build( mat );
}
CColMatrix &CColMatrix::operator=( const CColMatrix &mat )
{
build( mat );
return( *this );
}
void CColMatrix::build( const class CRowMatrix &mat )
{
_n = mat._n;
_m = mat._m;
_nz = mat._nz;
_asize = mat._nz;
reallocate();
int *c = new int[_m];
int i, j, ind, cl;
/* Count number of entries in each column to c. */
memset( c, 0, _m*sizeof(int) );
for( i = 0; i < _nz; i++ )
c[mat._col[i]]++;
/* Set up ptr. */
_ptr[0] = 0;
for( i = 1; i <= _m; i++ )
_ptr[i] = _ptr[i-1] + c[i-1];
/* Copy input to output using c as a row pointer. */
for( i = 0, j = 0; i < _n; i++ ) {
for( ; j < mat._ptr[i+1]; j++ ) {
cl = mat._col[j];
ind = _ptr[cl+1] - c[cl];
c[cl]--;
_row[ind] = i;
_val[ind] = mat._val[j];
}
}
delete [] c;
}
CColMatrix::CColMatrix( const class CRowMatrix &mat )
: _ptr(NULL), _row(NULL), _val(NULL)
{
build( mat );
}
CColMatrix &CColMatrix::operator=( const class CRowMatrix &mat )
{
build( mat );
return( *this );
}
void CColMatrix::build( const class CoordMatrix &mat )
{
_n = mat._n;
_m = mat._m;
_nz = mat._nz;
_asize = mat._nz;
reallocate();
int *c = new int[_m];
int i, start, cl;
/* Count number of entries in each column to c. */
memset( c, 0, _m*sizeof(int) );
for( i = 0; i < _nz; i++ )
c[mat._col[i]]++;
/* Set up ptr. */
_ptr[0] = 0;
for( i = 1; i <= _m; i++ )
_ptr[i] = _ptr[i-1] + c[i-1];
/* Copy input to output using c as a row pointer. */
for( i = 0; i < _nz; i++ ) {
cl = mat._col[i];
start = _ptr[cl];
c[cl]--;
_row[start+c[cl]] = mat._row[i];
_val[start+c[cl]] = mat._val[i];
}
delete [] c;
}
CColMatrix::CColMatrix( const class CoordMatrix &mat )
: _ptr(NULL), _row(NULL), _val(NULL)
{
build( mat );
}
CColMatrix &CColMatrix::operator=( const class CoordMatrix &mat )
{
build( mat );
return( *this );
}
CColMatrix::CColMatrix( const class Matrix &mat )
: _ptr(NULL), _row(NULL), _val(NULL)
{
const CRowMatrix *crmat;
const CColMatrix *ccmat;
const CoordMatrix *comat;
if( (crmat = dynamic_cast<const CRowMatrix *>(&mat)) != 0 )
build( *crmat );
else if( (ccmat = dynamic_cast<const CColMatrix *>(&mat)) != 0 )
build( *ccmat );
else if( (comat = dynamic_cast<const CoordMatrix *>(&mat)) != 0 )
build( *ccmat );
else
throw( ErrorUnimplemented( ERROR_LOCATION, "Couldn't convert unknown matrix type" ) );
}
CColMatrix &CColMatrix::operator=( const Matrix &mat )
{
const CRowMatrix *crmat;
const CColMatrix *ccmat;
const CoordMatrix *comat;
if( (crmat = dynamic_cast<const CRowMatrix *>(&mat)) != 0 )
build( *crmat );
else if( (ccmat = dynamic_cast<const CColMatrix *>(&mat)) != 0 )
build( *ccmat );
else if( (comat = dynamic_cast<const CoordMatrix *>(&mat)) != 0 )
build( *ccmat );
else
throw( ErrorUnimplemented( ERROR_LOCATION, "Couldn't convert unknown matrix type" ) );
return( *this );
}
CColMatrix::~CColMatrix()
{
free( _ptr );
free( _row );
free( _val );
}
void CColMatrix::resize( int n, int m )
{
free( _row );
free( _val );
_row = NULL;
_val = NULL;
if( _m != m ) {
int *tmp;
if( !(tmp = (int *)realloc( _ptr, (m+1)*sizeof(int) )) ) {
free( _ptr );
_ptr = NULL;
_n = _m = _nz = _asize = 0;
throw( ErrorNoMem( ERROR_LOCATION ) );
}
_ptr = tmp;
}
_n = n;
_m = m;
_nz = 0;
_asize = 0;
memset( _ptr, 0, (_m+1)*sizeof(int) );
}
void CColMatrix::clear( void )
{
_nz = 0;
_asize = 0;
memset( _ptr, 0, (_m+1)*sizeof(int) );
free( _row );
free( _val );
_row = NULL;
_val = NULL;
}
inline void CColMatrix::clear_no_check( int i, int j )
{
int a;
/* Search for the element */
for( a = _ptr[i]; a < _ptr[i+1]; a++ )
if( _row[a] == j )
break;
/* Do nothing if the element does not exist */
if( a == _ptr[i+1] )
return;
/* Move data */
int movesize = _nz-a-1;
memmove( &_val[a], &_val[a+1], movesize*sizeof(double) );
memmove( &_row[a], &_row[a+1], movesize*sizeof(int) );
/* Update pointers */
_nz--;
for( int a = i+1; a < _m+1; a++ )
_ptr[a]--;
}
void CColMatrix::clear_check( int i, int j )
{
if( i >= _n || j >= _m )
throw( ErrorRange( ERROR_LOCATION, i, _n, j, _m ) );
clear_no_check( i, j );
}
void CColMatrix::reserve( int size )
{
if( size > _asize ) {
_asize = size;
reallocate();
}
}
void CColMatrix::set_nz( int nz )
{
if( nz > _asize ) {
_asize = nz;
reallocate();
}
_nz = nz;
}
void CColMatrix::merge( CColMatrix &mat )
{
_n = mat._n;
_m = mat._m;
_nz = mat._nz;
_asize = mat._asize;
_ptr = mat._ptr;
_row = mat._row;
_val = mat._val;
mat._n = 0;
mat._m = 0;
mat._nz = 0;
mat._asize = 0;
mat.allocate();
}
void CColMatrix::order_ascending( void )
{
/* Sort each column. */
for( int i = 0; i < _m; i++ )
insertion_sort_iv( _row, _val, _ptr[i], _ptr[i+1] );
}
bool CColMatrix::check_ascending( void ) const
{
/* Check each column. */
for( int i = 0; i < _m; i++ ) {
int prev = 0;
for( int j = _ptr[i]; j < _ptr[i+1]; j++ ) {
if( j > _ptr[i] && prev > _row[j] )
return( false );
prev = _row[j];
}
}
return( true );
}
inline double CColMatrix::get_no_check( int i, int j ) const
{
for( int a = _ptr[j]; a < _ptr[j+1]; a++ )
if( _row[a] == i )
return( _val[a] );
return( 0.0 );
}
inline double &CColMatrix::set_no_check( int i, int j )
{
/* Use existing element if it exists */
for( int a = _ptr[j]; a < _ptr[j+1]; a++ )
if( _row[a] == i )
return( _val[a] );
/* Reserve new space if necessary */
if( _nz+1 > _asize )
reserve( _asize+_n );
/* Move existing data */
int movesize = _nz-_ptr[j+1];
memmove( &_val[_ptr[j+1]+1], &_val[_ptr[j+1]], movesize*sizeof(double) );
memmove( &_row[_ptr[j+1]+1], &_row[_ptr[j+1]], movesize*sizeof(int) );
/* Set new data */
_val[_ptr[j+1]] = 0.0;
_row[_ptr[j+1]] = i;
/* Update pointers */
_nz++;
for( int a = j+1; a < _m+1; a++ )
_ptr[a]++;
return( _val[_ptr[j+1]-1] );
}
double CColMatrix::get_check( int i, int j ) const
{
if( i >= _n || j >= _m )
throw( ErrorRange( ERROR_LOCATION, i, _n, j, _m ) );
return( get_no_check( i, j ) );
}
double &CColMatrix::set_check( int i, int j )
{
if( i >= _n || j >= _m )
throw( ErrorRange( ERROR_LOCATION, i, _n, j, _m ) );
return( set_no_check( i, j ) );
}
void CColMatrix::set_column( int j, int N, const int *row, const double *val )
{
if( j >= _m )
throw( ErrorRange( ERROR_LOCATION, j, _m ) );
else if( N > _n )
throw( Error( ERROR_LOCATION, "too many elements" ) );
/* Reserve new space if necessary */
int oldsize = _ptr[j+1] - _ptr[j];
if( _nz+N-oldsize > _asize )
reserve( _asize+_n );
/* Move existing data */
int offset = N-oldsize;
int movesize = _nz-_ptr[j+1];
memmove( &_val[_ptr[j+1]+offset], &_val[_ptr[j+1]], movesize*sizeof(double) );
memmove( &_row[_ptr[j+1]+offset], &_row[_ptr[j+1]], movesize*sizeof(int) );
/* Set new data */
int err = 0;
for( int a = 0; a < N; a++ ) {
for( int b = a+1; b < N; b++ ) {
if( row[a] == row[b] )
err = 2;
}
_val[_ptr[j]+a] = val[a];
if( (_row[_ptr[j]+a] = row[a]) >= _n )
err = 1;
}
/* Update pointers */
_nz += offset;
for( int a = j+1; a < _m+1; a++ )
_ptr[a] += offset;
if( err == 1 )
throw( Error( ERROR_LOCATION, "row index out of range" ) );
else if( err == 2 )
throw( Error( ERROR_LOCATION, "repeated row index" ) );
}
void CColMatrix::construct_add( int i, int j, double val )
{
/* Reserve new space if necessary */
if( _nz+1 > _asize )
reserve( _asize+_n );
/* Set new data */
_val[_ptr[j+1]] = val;
_row[_ptr[j+1]] = i;
/* Update the number of nonzeroes and only the next pointer */
_nz++;
_ptr[j+1]++;
if( j+1 < _m )
_ptr[j+2] = _ptr[j+1];
}
void CColMatrix::debug_print( std::ostream &os ) const
{
os << "n = " << _n << "\n";
os << "m = " << _m << "\n";
os << "nz = " << _nz << "\n";
os << "asize = " << _asize << "\n";
os << "ptr[] = {";
for( int i = 0; i < _m; i++ )
os << _ptr[i] << ", ";
os << _ptr[_m] << "}\n";
os << "row[] = {";
if( _nz <= 0 ) {
os << "}\n";
} else {
for( int i = 0; i < _nz-1; i++ )
os << _row[i] << ", ";
os << _row[_nz-1] << "}\n";
}
os << "val[] = {";
if( _nz <= 0 ) {
os << "}\n";
} else {
for( int i = 0; i < _nz-1; i++ )
os << _val[i] << ", ";
os << _val[_nz-1] << "}\n";
}
}
/* ************************************** *
* Matrix-Vector operations *
* ************************************** */
void CColMatrix::multiply_by_vector( Vector &x, const Vector &b ) const
{
if( b.size() != _m )
throw( ErrorDim( ERROR_LOCATION, "matrix dimension does not match vector" ) );
x.resize( _n );
x.clear();
double mult;
for( int i = 0; i < _n; i++ ) {
mult = b._val[i];
int end = _ptr[i+1];
for( int a = _ptr[i]; a < end; a++ )
x._val[_row[a]] += _val[a] * mult;
}
}
void CColMatrix::lower_unit_solve( Vector &x, const Vector &b ) const
{
// Make checks
if( _n != _m )
throw( ErrorDim( ERROR_LOCATION, "matrix not squrare" ) );
if( b.size() != _m )
throw( ErrorDim( ERROR_LOCATION, "matrix dimension does not match vector" ) );
x = b;
double mult;
for( int i = 0; i < _m; i++ ) {
mult = x[i];
for( int j = _ptr[i]; j < _ptr[i+1]; j++ )
x[_row[j]] -= _val[j] * mult;
}
}
void CColMatrix::upper_diag_solve( Vector &x, const Vector &b ) const
{
// Make checks
if( _n != _m )
throw( ErrorDim( ERROR_LOCATION, "matrix not square" ) );
if( b.size() != _m )
throw( ErrorDim( ERROR_LOCATION, "matrix dimension does not match vector" ) );
x = b;
double mult;
int32_t i, j;
for( i = _m-1; i >= 0; i-- ) {
j = _ptr[i+1]-1;
x[i] /= _val[j];
mult = x[i];
for( j--; j >= (int32_t)_ptr[i]; j-- )
x[_row[j]] -= _val[j] * mult;
}
}
