头文件：
/*
* Copyright (c) 2008-2011 Zhang Ming (M. Zhang), zmjerry@163.com
*
* This program 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 2 or any later version.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions are met:
*
* 1. Redistributions of source code must retain the above copyright notice,
* this list of conditions and the following disclaimer.
*
* 2. Redistributions in binary form must reproduce the above copyright
* notice, this list of conditions and the following disclaimer in the
* documentation and/or other materials provided with the distribution.
*
* This program 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. A copy of the GNU General Public License is available at:
* http://www.fsf.org/licensing/licenses
*/
/*****************************************************************************
* usingdeclare.h
*
* Usually used name declaration.
*
* Zhang Ming, 2010-08, Xi'an Jiaotong University.
*****************************************************************************/
#ifndef USINGDECLARE_H
#define USINGDECLARE_H
namespace splab
{
using std::cin;
using std::cout;
using std::cerr;
using std::endl;
using std::string;
using std::complex;
using std::ostream;
using std::istream;
using std::min;
using std::max;
using std::swap;
using std::ceil;
using std::abs;
using std::cos;
using std::sin;
using std::tan;
using std::acos;
using std::asin;
using std::atan;
using std::exp;
using std::log;
using std::log10;
using std::sqrt;
using std::pow;
}
// namespace splab
#endif
// USINGDECLARE_H/*
* Copyright (c) 2008-2011 Zhang Ming (M. Zhang), zmjerry@163.com
*
* This program 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 2 or any later version.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions are met:
*
* 1. Redistributions of source code must retain the above copyright notice,
* this list of conditions and the following disclaimer.
*
* 2. Redistributions in binary form must reproduce the above copyright
* notice, this list of conditions and the following disclaimer in the
* documentation and/or other materials provided with the distribution.
*
* This program 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. A copy of the GNU General Public License is available at:
* http://www.fsf.org/licensing/licenses
*/
/*****************************************************************************
* constants.h
*
* Some constants often used in numeric computing.
*
* Zhang Ming, 2010-01, Xi'an Jiaotong University.
*****************************************************************************/
#ifndef CONSTANTS_H
#define CONSTANTS_H
namespace splab
{
const double EPS = 2.220446049250313e-016;
const double PI = 3.141592653589793;
const double TWOPI = 6.283185307179586;
const double HALFPI = 1.570796326794897;
const double D2R = 0.017453292519943;
const double EXP = 2.718281828459046;
const double RT2 = 1.414213562373095;
const double IRT2 = 0.707106781186547;
const int FORWARD = 1;
const int INVERSE = 0;
const int MAXTERM = 20;
const int INITSIZE = 20;
const int EXTFACTOR = 2;
}
// namespace splab
#endif
// CONSTANTS_H/*
* Copyright (c) 2008-2011 Zhang Ming (M. Zhang), zmjerry@163.com
*
* This program 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 2 or any later version.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions are met:
*
* 1. Redistributions of source code must retain the above copyright notice,
* this list of conditions and the following disclaimer.
*
* 2. Redistributions in binary form must reproduce the above copyright
* notice, this list of conditions and the following disclaimer in the
* documentation and/or other materials provided with the distribution.
*
* This program 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. A copy of the GNU General Public License is available at:
* http://www.fsf.org/licensing/licenses
*/
/*****************************************************************************
* vector.h
*
* Class template of vector which is designed for basic linear algebra
* operations such as:
* v + k k + v v += k v1 + v2 v1 += v2
* v - k k - v v -= k v1 - v2 v1 -= v2
* v * k k * v v *= k v1 * v2 v1 *= v2
* v / k k / v v /= k v1 / v2 v1 /= v2
* mum, min, max swap reverse
* norm dotProd
* These operators and functions can be applied to both real vector and
* complex vector.
*
* The class also provides the basic math functions such as:
* cos sin tan acos asin atan
* abs exp log log10 sqrt pow
* This should include "matrixmath.h" file.
*
* When debugging, use #define BOUNDS_CHECK above your "#include vector.h"
* line. When done debugging, comment out #define BOUNDS_CHECK for better
* performance.
*
* Zhang Ming, 2010-01 (revised 2010-12), Xi'an Jiaotong University.
*****************************************************************************/
#ifndef VECTOR_H
#define VECTOR_H
#include <iostream>
#include <cassert>
#include <cmath>
#include <complex>
#include <usingdeclare.h>
#include <constants.h>
namespace splab
{
template <typename Type>
class Vector
{
public:
typedef Type* iterator;
typedef const Type* const_iterator;
// constructors and destructor
Vector();
Vector( const Vector<Type> &v );
Vector( int length, const Type &x = Type(0) );
Vector( int length, const Type *array );
~Vector();
// assignments
Vector<Type>& operator=( const Vector<Type> &v );
Vector<Type>& operator=( const Type &x );
// accessors
Type& operator[]( int i );
const Type& operator[]( int i ) const;
Type& operator()( int i );
const Type& operator()( int i ) const;
// iterators
iterator begin();
const_iterator begin() const;
iterator end();
const_iterator end() const;
// type conversion
operator Type*();
operator const Type*() const;
// others
int size() const;
int dim() const;
Vector<Type>& resize( int length );
// computed assignment
Vector<Type>& operator+=( const Type& );
Vector<Type>& operator-=( const Type& );
Vector<Type>& operator*=( const Type& );
Vector<Type>& operator/=( const Type& );
Vector<Type>& operator+=( const Vector<Type>& );
Vector<Type>& operator-=( const Vector<Type>& );
Vector<Type>& operator*=( const Vector<Type>& );
Vector<Type>& operator/=( const Vector<Type>& );
private:
// data pointer for 0-offset indexing
Type *pv0;
// data pointer for 1-offset indexing
Type *pv1;
// the row number of vector
int nRow;
void init( int length );
void copyFromArray( const Type *v );
void setByScalar( const Type &x );
void destroy();
};
// class Vector
// input and output
template<typename Type>
ostream& operator<<( ostream&, const Vector<Type>& );
template<typename Type>
istream& operator>>( istream&, Vector<Type>& );
// arithmetic operators
template<typename Type>
Vector<Type> operator-( const Vector<Type>& );
template<typename Type>
Vector<Type> operator+( const Vector<Type>&, const Type& );
template<typename Type>
Vector<Type> operator+( const Type&, const Vector<Type>& );
template<typename Type>
Vector<Type> operator+( const Vector<Type>&, const Vector<Type>& );
template<typename Type>
Vector<Type> operator-( const Vector<Type>&, const Type& );
template<typename Type>
Vector<Type> operator-( const Type&, const Vector<Type>& );
template<typename Type>
Vector<Type> operator-( const Vector<Type>&, const Vector<Type>& );
template<typename Type>
Vector<Type> operator*( const Vector<Type>&, const Type& );
template<typename Type>
Vector<Type> operator*( const Type&, const Vector<Type>& );
template<typename Type>
Vector<Type> operator*( const Vector<Type>&, const Vector<Type>& );
template<typename Type>
Vector<Type> operator/( const Vector<Type>&, const Type& );
template<typename Type>
Vector<Type> operator/( const Type&, const Vector<Type>& );
template<typename Type>
Vector<Type> operator/( const Vector<Type>&, const Vector<Type>& );
// dot product
template<typename Type>
Type dotProd( const Vector<Type>&, const Vector<Type>& );
template<typename Type> complex<Type>
dotProd( const Vector<complex<Type> >&, const Vector<complex<Type> >& );
// utilities
template<typename Type> Type sum( const Vector<Type>& );
template<typename Type> Type min( const Vector<Type>& );
template<typename Type> Type max( const Vector<Type>& );
template<typename Type> Type norm( const Vector<Type>& );
template<typename Type> Type norm( const Vector<complex<Type> >& );
template<typename Type> void swap( Vector<Type>&, Vector<Type>& );
template<typename Type> Vector<Type> linspace( Type, Type, int );
template<typename Type> Vector<Type> abs( const Vector<complex<Type> >& );
template<typename Type> Vector<Type> arg( const Vector<complex<Type> >& );
template<typename Type> Vector<Type> real( const Vector<complex<Type> >& );
template<typename Type> Vector<Type> imag( const Vector<complex<Type> >& );
template<typename Type>
Vector<complex<Type> > complexVector( const Vector<Type>& );
template<typename Type>
Vector<complex<Type> > complexVector( const Vector<Type>&,
const Vector<Type>& );
#include <vector-impl.h>
}
// namespace splab
#endif
// VECTOR_H实现文件：
/*
* Copyright (c) 2008-2011 Zhang Ming (M. Zhang), zmjerry@163.com
*
* This program 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 2 or any later version.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions are met:
*
* 1. Redistributions of source code must retain the above copyright notice,
* this list of conditions and the following disclaimer.
*
* 2. Redistributions in binary form must reproduce the above copyright
* notice, this list of conditions and the following disclaimer in the
* documentation and/or other materials provided with the distribution.
*
* This program 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. A copy of the GNU General Public License is available at:
* http://www.fsf.org/licensing/licenses
*/
/*****************************************************************************
* vector-impl.h
*
* Implementation for Vector class.
*
* Zhang Ming, 2010-01 (revised 2010-12), Xi'an Jiaotong University.
*****************************************************************************/
/**
* initialize
*/
template <typename Type>
void Vector<Type>::init( int length )
{
assert( pv0 == NULL );
pv0 = new Type[length];
assert( pv0 != NULL );
pv1 = pv0 - 1;
nRow = length;
}
/**
* copy vector from normal array
*/
template <typename Type>
inline void Vector<Type>::copyFromArray( const Type *v )
{
for( int i=0; i<nRow; ++i )
pv0[i] = v[i];
}
/**
* set vector by a scalar
*/
template <typename Type>
inline void Vector<Type>::setByScalar( const Type &x )
{
for( int i=0; i<nRow; ++i )
pv0[i] = x;
}
/**
* destroy the vector
*/
template <typename Type>
void Vector<Type>::destroy()
{
if( pv0 == NULL )
return;
delete []pv0;
pv0 = NULL;
pv1 = NULL;
}
/**
* constructors and destructor
*/
template <typename Type>
Vector<Type>::Vector()
: pv0(0), pv1(0), nRow(0)
{
}
template <typename Type>
Vector<Type>::Vector( const Vector<Type> &v )
: pv0(0), pv1(0), nRow(0)
{
init( v.nRow );
copyFromArray( v.pv0 );
}
template <typename Type>
Vector<Type>::Vector( int length, const Type &x )
: pv0(0), pv1(0), nRow(0)
{
init( length );
setByScalar( x );
}
template <typename Type>
Vector<Type>::Vector( int length, const Type *array )
: pv0(0), pv1(0), nRow(0)
{
init( length );
copyFromArray( array );
}
template <typename Type>
Vector<Type>::~Vector()
{
destroy();
}
/**
* overload evaluate operator= from vector to vector
*/
template <typename Type>
Vector<Type>& Vector<Type>::operator=( const Vector<Type> &v )
{
if( pv0 == v.pv0 )
return *this;
if( nRow == v.nRow )
copyFromArray( v.pv0 );
else
{
destroy();
init( v.nRow );
copyFromArray( v.pv0 );
}
return *this;
}
/**
* overload evaluate operator= from scalar to vector
*/
template <typename Type>
inline Vector<Type>& Vector<Type>::operator=( const Type &x )
{
setByScalar( x );
return *this;
}
/**
* overload operator [] for 0-offset access
*/
template <typename Type>
inline Type& Vector<Type>::operator[]( int i )
{
#ifdef BOUNDS_CHECK
assert( 0 <= i );
assert( i < nRow );
#endif
return pv0[i];
}
template <typename Type>
inline const Type& Vector<Type>::operator[]( int i ) const
{
#ifdef BOUNDS_CHECK
assert( 0 <= i );
assert( i < nRow );
#endif
return pv0[i];
}
/**
* overload operator () for 1-offset access
*/
template <typename Type>
inline Type& Vector<Type>::operator()( int i )
{
#ifdef BOUNDS_CHECK
assert( 1 <= i );
assert( i <= nRow );
#endif
return pv1[i];
}
template <typename Type>
inline const Type& Vector<Type>::operator()( int i ) const
{
#ifdef BOUNDS_CHECK
assert( 1 <= i );
assert( i <= nRow );
#endif
return pv1[i];
}
/**
* iterators
*/
template <typename Type>
inline typename Vector<Type>::iterator Vector<Type>::begin()
{
return pv0;
}
template <typename Type>
inline typename Vector<Type>::const_iterator Vector<Type>::begin() const
{
return pv0;
}
template <typename Type>
inline typename Vector<Type>::iterator Vector<Type>::end()
{
return pv0 + nRow;
}
template <typename Type>
inline typename Vector<Type>::const_iterator Vector<Type>::end() const
{
return pv0 + nRow;
}
/**
* type conversion functions
*/
template <typename Type>
inline Vector<Type>::operator Type*()
{
return pv0;
}
template <typename Type>
inline Vector<Type>::operator const Type*() const
{
return pv0;
}
/**
* get the vector's total size
*/
template <typename Type>
inline int Vector<Type>::size() const
{
return nRow;
}
/**
* get the vector's dimension
*/
template <typename Type>
inline int Vector<Type>::dim() const
{
return nRow;
}
/**
* reallocate vector's size
*/
template <typename Type>
Vector<Type>& Vector<Type>::resize( int length )
{
if( nRow == length )
return *this;
destroy();
init( length );
return *this;
}
/**
* compound assignment operators +=
*/
template <typename Type>
Vector<Type>& Vector<Type>::operator+=( const Type &x )
{
iterator itr = (*this).begin();
while( itr != (*this).end() )
*itr++ += x;
return *this;
}
template <typename Type>
Vector<Type>& Vector<Type>::operator+=( const Vector<Type> &rhs )
{
assert( nRow == rhs.dim() );
iterator itrL = (*this).begin();
const_iterator itrR = rhs.begin();
while( itrL != (*this).end() )
*itrL++ += *itrR++;
return *this;
}
/**
* compound assignment operators -=
*/
template <typename Type>
Vector<Type>& Vector<Type>::operator-=( const Type &x )
{
iterator itr = (*this).begin();
while( itr != (*this).end() )
*itr++ -= x;
return *this;
}
template <typename Type>
Vector<Type>& Vector<Type>::operator-=( const Vector<Type> &rhs )
{
assert( nRow == rhs.dim() );
iterator itrL = (*this).begin();
const_iterator itrR = rhs.begin();
while( itrL != (*this).end() )
*itrL++ -= *itrR++;
return *this;
}
/**
* compound assignment operators *=
*/
template <typename Type>
Vector<Type>& Vector<Type>::operator*=( const Type &x )
{
iterator itr = (*this).begin();
while( itr != (*this).end() )
*itr++ *= x;
return *this;
}
template <typename Type>
Vector<Type>& Vector<Type>::operator*=( const Vector<Type> &rhs )
{
assert( nRow == rhs.dim() );
iterator itrL = (*this).begin();
const_iterator itrR = rhs.begin();
while( itrL != (*this).end() )
*itrL++ *= *itrR++;
return *this;
}
/**
* compound assignment operators /=
*/
template <typename Type>
Vector<Type>& Vector<Type>::operator/=( const Type &x )
{
iterator itr = (*this).begin();
while( itr != (*this).end() )
*itr++ /= x;
return *this;
}
template <typename Type>
Vector<Type>& Vector<Type>::operator/=( const Vector<Type> &rhs )
{
assert( nRow == rhs.dim() );
iterator itrL = (*this).begin();
const_iterator itrR = rhs.begin();
while( itrL != (*this).end() )
*itrL++ /= *itrR++;
return *this;
}
/**
* Overload the output stream function.
*/
template <typename Type>
ostream& operator<<( ostream &out, const Vector<Type> &v )
{
int N = v.dim();
out << "size: " << N << " by 1" << "\n";
for( int i=0; i<N; ++i )
out << v[i] << " " << "\n";
return out;
}
/**
* Overload the input stream function.
*/
template <typename Type>
istream& operator>>( istream &in, Vector<Type> &v )
{
int N;
in >> N;
if( !( N == v.dim() ) )
v.resize( N );
for( int i=0; i<N; ++i )
in >> v[i];
return in;
}
/**
* get negative vector
*/
template <typename Type>
Vector<Type> operator-( const Vector<Type> &v )
{
Vector<Type> tmp( v.dim() );
typename Vector<Type>::iterator itrL = tmp.begin();
typename Vector<Type>::const_iterator itrR = v.begin();
while( itrL != tmp.end() )
*itrL++ = -(*itrR++);
return tmp;
}
/**
* vector-scalar addition.
*/
template <typename Type>
inline Vector<Type> operator+( const Vector<Type> &v, const Type &x )
{
Vector<Type> tmp( v );
return tmp += x;
}
template <typename Type>
inline Vector<Type> operator+( const Type &x, const Vector<Type> &v )
{
return v+x;
}
/**
* vector-scalar substraction.
*/
template <typename Type>
inline Vector<Type> operator-( const Vector<Type> &v, const Type &x )
{
Vector<Type> tmp( v );
return tmp -= x;
}
template <typename Type>
inline Vector<Type> operator-( const Type &x, const Vector<Type> &v )
{
Vector<Type> tmp( v );
return -tmp += x;
}
/**
* vector-scalar multiplication.
*/
template <typename Type>
inline Vector<Type> operator*( const Vector<Type> &v, const Type &x )
{
Vector<Type> tmp( v );
return tmp *= x;
}
template <typename Type>
inline Vector<Type> operator*( const Type &x, const Vector<Type> &v )
{
return v*x;
}
/**
* vector-scalar division.
*/
template <typename Type>
inline Vector<Type> operator/( const Vector<Type> &v, const Type &x )
{
Vector<Type> tmp( v );
return tmp /= x;
}
template <typename Type>
inline Vector<Type> operator/( const Type &x, const Vector<Type> &v )
{
int N = v.dim();
Vector<Type> tmp( N );
for( int i=0; i<N; ++i )
tmp[i] = x / v[i];
return tmp;
}
/**
* vector-vector addition.
*/
template <typename Type>
inline Vector<Type> operator+( const Vector<Type> &v1, const Vector<Type> &v2 )
{
Vector<Type> tmp( v1 );
return tmp += v2;
}
/**
* vector-vector substraction.
*/
template <typename Type>
inline Vector<Type> operator-( const Vector<Type> &v1, const Vector<Type> &v2 )
{
Vector<Type> tmp( v1 );
return tmp -= v2;
}
/**
* vector-vector multiplication.
*/
template <typename Type>
inline Vector<Type> operator*( const Vector<Type> &v1, const Vector<Type> &v2 )
{
Vector<Type> tmp( v1 );
return tmp *= v2;
}
/**
* vector-vector division.
*/
template <typename Type>
inline Vector<Type> operator/( const Vector<Type> &v1, const Vector<Type> &v2 )
{
Vector<Type> tmp( v1 );
return tmp /= v2;
}
/**
* Inner product for vectors.
*/
template <typename Type>
Type dotProd( const Vector<Type> &v1, const Vector<Type> &v2 )
{
assert( v1.dim() == v2.dim() );
Type sum = 0;
typename Vector<Type>::const_iterator itr1 = v1.begin();
typename Vector<Type>::const_iterator itr2 = v2.begin();
while( itr1 != v1.end() )
sum += (*itr1++) * (*itr2++);
return sum;
}
/**
* Inner product for vectors.
*/
template <typename Type>
complex<Type> dotProd( const Vector<complex<Type> > &v1,
const Vector<complex<Type> > &v2 )
{
assert( v1.dim() == v2.dim() );
complex<Type> sum = 0;
typename Vector<complex<Type> >::const_iterator itr1 = v1.begin();
typename Vector<complex<Type> >::const_iterator itr2 = v2.begin();
while( itr1 != v1.end() )
sum += (*itr1++) * conj(*itr2++);
return sum;
}
/**
* Vector's sum.
*/
template <typename Type>
Type sum( const Vector<Type> &v )
{
Type sum = 0;
typename Vector<Type>::const_iterator itr = v.begin();
while( itr != v.end() )
sum += *itr++;
return sum;
}
/**
* Minimum value of vector.
*/
template <typename Type>
Type min( const Vector<Type> &v )
{
Type m = v[0];
for( int i=1; i<v.size(); ++i )
if( m > v[i] )
m = v[i];
return m;
}
/**
* Maximum value of vector.
*/
template <typename Type>
Type max( const Vector<Type> &v )
{
Type M = v[0];
for( int i=1; i<v.size(); ++i )
if( M < v[i] )
M = v[i];
return M;
}
/**
* Vector's norm in Euclidean space.
*/
template <typename Type>
Type norm( const Vector<Type> &v )
{
Type sum = 0;
typename Vector<Type>::const_iterator itr = v.begin();
while( itr != v.end() )
{
sum += (*itr) * (*itr);
itr++;
}
return Type(sqrt(1.0*sum));
}
/**
* Vector's norm in Euclidean space.
*/
template <typename Type>
Type norm( const Vector<complex<Type> > &v )
{
Type sum = 0;
typename Vector<complex<Type> >::const_iterator itr = v.begin();
while( itr != v.end() )
sum += norm(*itr++);
return Type(sqrt(1.0*sum));
}
/**
* return vector's reversion
*/
template <typename Type>
void swap( Vector<Type> &lhs, Vector<Type> &rhs )
{
typename Vector<Type>::iterator itrL = lhs.begin(),
itrR = rhs.begin();
while( itrL != lhs.end() )
std::swap( *itrL++, *itrR++ );
}
/**
* Generates a vector of n points linearly spaced between and
* including a and b.
*/
template <typename Type>
Vector<Type> linspace( Type a, Type b, int n )
{
if( n < 1 )
return Vector<Type>();
else if( n == 1 )
return Vector<Type>( 1, a );
else
{
Type dx = (b-a) / (n-1);
Vector<Type> tmp(n);
for( int i=0; i<n; ++i )
tmp[i] = a + i*dx;
return tmp;
}
}
/**
* Get magnitude of a complex vector.
*/
template <typename Type>
Vector<Type> abs( const Vector< complex<Type> > &v )
{
Vector<Type> tmp( v.dim() );
typename Vector<Type>::iterator itrL = tmp.begin();
typename Vector< complex<Type> >::const_iterator itrR = v.begin();
while( itrL != tmp.end() )
*itrL++ = abs(*itrR++);
return tmp;
}
/**
* Get angle of a complex vector.
*/
template <typename Type>
Vector<Type> arg( const Vector< complex<Type> > &v )
{
Vector<Type> tmp( v.dim() );
typename Vector<Type>::iterator itrL = tmp.begin();
typename Vector< complex<Type> >::const_iterator itrR = v.begin();
while( itrL != tmp.end() )
*itrL++ = arg(*itrR++);
return tmp;
}
/**
* Get real part of a complex vector.
*/
template <typename Type>
Vector<Type> real( const Vector< complex<Type> > &v )
{
Vector<Type> tmp( v.dim() );
typename Vector<Type>::iterator itrL = tmp.begin();
typename Vector< complex<Type> >::const_iterator itrR = v.begin();
while( itrL != tmp.end() )
*itrL++ = (*itrR++).real();
return tmp;
}
/**
* Get imaginary part of a complex vector.
*/
template <typename Type>
Vector<Type> imag( const Vector< complex<Type> > &v )
{
Vector<Type> tmp( v.dim() );
typename Vector<Type>::iterator itrL = tmp.begin();
typename Vector< complex<Type> >::const_iterator itrR = v.begin();
while( itrL != tmp.end() )
*itrL++ = (*itrR++).imag();
return tmp;
}
/**
* Convert real vector to complex vector.
*/
template <typename Type>
Vector<complex<Type> > complexVector( const Vector<Type> &rv )
{
int N = rv.dim();
Vector<complex<Type> > cv( N );
typename Vector<complex<Type> >::iterator itrL = cv.begin();
typename Vector<Type>::const_iterator itrR = rv.begin();
while( itrR != rv.end() )
*itrL++ = *itrR++;
return cv;
}
template <typename Type>
Vector<complex<Type> > complexVector( const Vector<Type> &vR,
const Vector<Type> &vI )
{
int N = vR.dim();
assert( N == vI.dim() );
Vector<complex<Type> > cv( N );
typename Vector<complex<Type> >::iterator itrC = cv.begin();
typename Vector<Type>::const_iterator itrR = vR.begin(),
itrI = vI.begin();
while( itrC != cv.end() )
*itrC++ = complex<Type>( *itrR++, *itrI++ );
return cv;
}测试代码：
/*****************************************************************************
* vector_test.cpp
*
* Vector class testing.
*
* Zhang Ming, 2010-01 (revised 2010-12), Xi'an Jiaotong University.
*****************************************************************************/
#define BOUNDS_CHECK
#include <iostream>
#include <iomanip>
#include <complex>
#include <vector.h>
using namespace std;
using namespace splab;
const int M = 3;
void display( const int *p, int length )
{
for( int i=0; i<length; ++i )
cout << p[i] << "\t" ;
cout << endl;
}
int main()
{
int k;
int arrays[3] = {1,2,3};
Vector<int> v1( M,arrays );
k = 1;
Vector<int> v2( M,k );
Vector<int> v3( M );
k = 0;
v3 = k;
Vector<int> v4( v1 );
cout << "vector v1 : " << v1 << endl;
cout << "vector v2 : " << v2 << endl;
cout << "vector v3 : " << v3 << endl;
cout << "vector v4 : " << v4 << endl;
display( v4, M );
cout << endl;
v4.resize( 5 );
Vector<int>::iterator itr = v4.begin();
while( itr != v4.end() )
*itr++ = 1;
cout << "new vector v4 : " << v4 << endl;
v4 = v1;
cout << "new vector v4 : " << v4 << endl;
k = 2;
v3 = k+v1;
cout << "v3 = k + v1 : " << v3 << endl;
v3 += k;
cout << "v3 += k : " << v3 << endl;
v3 = v1-k;
cout << "v3 = v1 - k : " << v3 << endl;
v3 = k-v1;
cout << "v3 = k - v1 : " << v3 << endl;
v3 -= k;
cout << "v3 -= k : " << v3 << endl;
v3 = k*v1;
cout << "v3 = k * v1 : " << v3 << endl;
v3 *= k;
cout << "v3 *= k : " << v3 << endl;
v3 = v1/k;
cout << "v3 = v1 / k : " << v3 << endl;
v3 = k/v1;
cout << "v3 = k / v1 : " << v3 << endl;
v3 /= k;
cout << "v3 /= k : " << v3 << endl;
v3 = v1+v2;
cout << "v3 = v1 + v2 : " << v3 << endl;
v3 += v1;
cout << "v3 += v1 : " << v3 << endl;
v3 = v1-v2;
cout << "v3 = v1 - v2 : " << v3 << endl;
v3 -= v1;
cout << "v3 -= v1 : " << v3 << endl;
v3 = v1*v2;
cout << "v3 = v1 * v2 : " << v3 << endl;
v3 *= v1;
cout << "v3 *= v1 : " << v3 << endl;
v3 = v1/v2;
cout << "v3 = v1 / v2 : " << v3 << endl;
v3 /= v1;
cout << "v3 /= v1 : " << v3 << endl;
cout << "minimum element of v1 : " << min(v1) << endl << endl;
cout << "maximum element of v1 : " << max(v1) << endl << endl;
cout << "L2 norm of v3 : " << norm( v3 ) << endl << endl;
cout << "inner product of v1 and v2 : " << dotProd( v1, v2 ) << endl << endl;
complex<double> z = -1.0;
Vector< complex<double> > v( M );
v[0] = polar( 1.0,PI/4 );
v[1] = polar( 1.0,PI );
v[2] = complex<double>( 1.0,-1.0 );
Vector< complex<double> > u = v*z;
cout << setiosflags(ios::fixed) << setprecision(4);
cout << "convert from real to complex vector: " << complexVector(v3) << endl;
cout << "convert from real to complex vector: " << complexVector(v3,-v1) << endl;
cout << "complex vector v : " << v << endl;
cout << "complex vector u = -v : " << u << endl;
cout << "norm of coplex vector v : " << norm(v) << endl << endl;
cout << "dot product of complex vector v and 1+u: " << dotProd(v,u-z) << endl << endl;
int N = 5;
Vector<double> x = linspace( 0.0, TWOPI, N );
Vector< complex<float> > cv(N);
for( int i=0; i<N; ++i )
cv[i] = complex<float>( float(sin(x[i])), float(cos(x[i])) );
cout << "Complex vector vc : " << cv << endl;
cout << "Absolute of vc : " << abs(cv) << endl;
cout << "Angle of vc : " << arg(cv) << endl;
cout << "Real part of vc : " << real(cv) << endl;
cout << "Imaginary part of vc : " << imag(cv) << endl;
cout << resetiosflags(ios::fixed);
Vector< Vector<double> > v2d1( M );
for( int i=0; i<M; ++i )
{
v2d1[i].resize( M );
for( int j=0; j<M; ++j )
v2d1[i][j] = double( i+j );
}
cout << "two dimension vector v2d1 : " << endl;
Vector< Vector<double> >::const_iterator itrD2 = v2d1.begin();
int rowNum = 0;
while( itrD2 != v2d1.end() )
cout << "the " << rowNum++ << "th row : " << *itrD2++ << endl;
Vector< Vector<double> > v2d2 = v2d1+v2d1;
cout << "two dimension vector v2d2 = v2d1 + v2d1 : " << v2d2;
return 0;
}运行结果：
size: 5 by 1
0.0000
1.5708
3.1416
4.7124
6.2832
sin of x : size: 5 by 1
0.0000
1.0000
0.0000
-1.0000
-0.0000
cos of x : size: 5 by 1
1.0000
0.0000
-1.0000
-0.0000
1.0000
tan of x : size: 5 by 1
0.0000
16331778728383844.0000
-0.0000
5443926242794615.0000
-0.0000
asin of x : size: 5 by 1
0.0000
nan
nan
nan
nan
acos of x : size: 5 by 1
1.5708
nan
nan
nan
nan
atan of x : size: 5 by 1
0.0000
1.0039
1.2626
1.3617
1.4130
exp of x : size: 5 by 1
1.0000
4.8105
23.1407
111.3178
535.4917
log of x : size: 5 by 1
-inf
0.4516
1.1447
1.5502
1.8379
log10 of x : size: 5 by 1
-inf
0.1961
0.4971
0.6732
0.7982
sqrt of x : size: 5 by 1
0.0000
1.2533
1.7725
2.1708
2.5066
pow of x : size: 5 by 1
1.0000
2.0327
36.4622
1487.9012
103540.9204
pow of x : size: 5 by 1
0.0000
2.4674
9.8696
22.2066
39.4784
pow of x : size: 5 by 1
1.0000
2.9707
8.8250
26.2162
77.8802
The standard normal distribution : size: 5 by 1
0.0604
0.0633
0.0625
0.0578
0.0503
Process returned 0 (0x0) execution time : 0.078 s
Press any key to continue.
1777
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
