头文件:
/*
* 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
*/
/*****************************************************************************
* matrix.h
*
* Class template of matrix which is designed for basic linear algebra
* operations such as:
* A + x x + A A += x A1 + A2 A1 += A2
* A - x x - A A -= x A1 - A2 A1 -= A2
* A * x x * A A *= x A1 * A2
* A / x x / A A /= x
* mum, min, max mean swap
* eye, diag, trT trH norm
* trMult multTr elemMult elemMultEq
* elemDivd elemDivdEq
* These operators and functions can be applied to both real matrix and
* complex matrix.
*
* 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 matrix.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 MATRIX_H
#define MATRIX_H
#include <vector.h>
namespace splab
{
template <typename Type>
class Matrix
{
public:
// constructors and destructor
Matrix();
Matrix( const Matrix<Type> &A );
Matrix( int rows, int columns, const Type &x = Type(0) );
Matrix( int rows, int columns, const Type *v );
~Matrix();
// assignments
Matrix<Type>& operator=( const Matrix<Type> &A );
Matrix<Type>& operator=( const Type &x );
// accessors
Type* operator[]( int i );
const Type* operator[]( int i ) const;
Type& operator()( int row, int column );
const Type& operator()( int row, int column ) const;
// type conversion
operator Type*();
operator const Type*() const;
// operator Type**();
// operator const Type**() const;
// others
long size() const;
int dim( int dimension ) const;
int rows() const;
int cols() const;
Matrix<Type>& resize( int rows, int columns );
Vector<Type> getRow( int row ) const;
Vector<Type> getColumn( int column ) const;
void setRow( const Vector<Type> &v, int row );
void setColumn( const Vector<Type> &v, int column );
// computed assignment
Matrix<Type>& operator+=( const Type& );
Matrix<Type>& operator+=( const Matrix<Type>& );
Matrix<Type>& operator-=( const Type& );
Matrix<Type>& operator-=( const Matrix<Type>& );
Matrix<Type>& operator*=( const Type& );
// WARNING: element-by-element
Matrix<Type>& operator*=( const Matrix<Type>& );
Matrix<Type>& operator/=( const Type& );
// WARNING: element-by-element
Matrix<Type>& operator/=( const Matrix<Type>& );
private:
// 0-based and 1-based data pointer
Type *pv0, *pv1;
// 0-based and 1-based row pointer's pointer
Type **prow0, **prow1;
// row number, column number and total number
int nRow;
int nColumn;
long nTotal;
void init( int rows, int columns );
void copyFromArray( const Type *v );
void setByScalar( const Type &x );
void destroy();
};
// class Matrix
// input and output
template<typename Type>
ostream& operator<<( ostream&, const Matrix<Type>& );
template<typename Type>
istream& operator>>( istream&, Matrix<Type>& );
// arithmetic operators
template<typename Type>
Matrix<Type> operator-( const Matrix<Type>& );
template<typename Type>
Matrix<Type> operator+( const Matrix<Type>&, const Type& );
template<typename Type>
Matrix<Type> operator+( const Type&, const Matrix<Type>& );
template<typename Type>
Matrix<Type> operator+( const Matrix<Type>&, const Matrix<Type>& );
template<typename Type>
Matrix<Type> operator-( const Matrix<Type>&, const Type& );
template<typename Type>
Matrix<Type> operator-( const Type&, const Matrix<Type>& );
template<typename Type>
Matrix<Type> operator-( const Matrix<Type>&, const Matrix<Type>& );
template<typename Type>
Matrix<Type> operator*( const Matrix<Type>&, const Type& );
template<typename Type>
Matrix<Type> operator*( const Type&, const Matrix<Type>& );
template<typename Type>
Matrix<Type> operator*( const Matrix<Type>&, const Matrix<Type>& );
template<typename Type>
Vector<Type> operator*( const Matrix<Type>&, const Vector<Type>& );
template<typename Type>
Matrix<Type> operator/( const Matrix<Type>&, const Type& );
template<typename Type>
Matrix<Type> operator/( const Type&, const Matrix<Type>& );
// optimizied version of multiplication
template<typename Type>
Matrix<Type>& optMult( const Matrix<Type>&, const Matrix<Type>&,
Matrix<Type>& );
template<typename Type>
Vector<Type>& optMult( const Matrix<Type>&, const Vector<Type>&,
Vector<Type>& );
// element-by-element multiplication and division
template<typename Type>
Matrix<Type> elemMult( const Matrix<Type>&, const Matrix<Type>& );
template<typename Type>
Matrix<Type> elemDivd( const Matrix<Type>&, const Matrix<Type>& );
template<typename Type>
Matrix<Type>& elemMultEq( Matrix<Type>&, const Matrix<Type>& );
template<typename Type>
Matrix<Type>& elemDivdEq( Matrix<Type>&, const Matrix<Type>& );
// transpose and conjugate transpose
template<typename Type> Matrix<Type> trT( const Matrix<Type>& );
template<typename Type> Matrix<Type> trH( const Matrix<Type>& );
// transpose multiplication
template<typename Type>
Matrix<Type> trMult( const Matrix<Type>&, const Matrix<Type>& );
template<typename Type>
Vector<Type> trMult( const Matrix<Type>&, const Vector<Type>& );
template<typename Type>
Matrix<Type> multTr( const Matrix<Type>&, const Matrix<Type>& );
template<typename Type>
Matrix<Type> multTr( const Vector<Type>&, const Vector<Type>& );
template<typename Type>
Matrix<complex<Type> > trMult( const Matrix<complex<Type> >&,
const Matrix<complex<Type> >& );
template<typename Type>
Vector<complex<Type> > trMult( const Matrix<complex<Type> >&,
const Vector<complex<Type> >& );
template<typename Type>
Matrix<complex<Type> > multTr( const Matrix<complex<Type> >&,
const Matrix<complex<Type> >& );
template<typename Type>
Matrix<complex<Type> > multTr( const Vector<complex<Type> >&,
const Vector<complex<Type> >& );
// unit and diagonal matrix
template<typename Type> Matrix<Type> eye( int, const Type& );
template<typename Type> Vector<Type> diag( const Matrix<Type>& );
template<typename Type> Matrix<Type> diag( const Vector<Type>& );
// utilities
template<typename Type> Type norm( const Matrix<Type>& );
template<typename Type> Type norm( const Matrix<complex<Type> >& );
template<typename Type> void swap( Matrix<Type>&, Matrix<Type>& );
template<typename Type> Vector<Type> sum( const Matrix<Type>& );
template<typename Type> Vector<Type> min( const Matrix<Type>& );
template<typename Type> Vector<Type> max( const Matrix<Type>& );
template<typename Type> Vector<Type> mean( const Matrix<Type>& );
template<typename Type> Matrix<Type> abs( const Matrix<complex<Type> >& );
template<typename Type> Matrix<Type> arg( const Matrix<complex<Type> >& );
template<typename Type> Matrix<Type> real( const Matrix<complex<Type> >& );
template<typename Type> Matrix<Type> imag( const Matrix<complex<Type> >& );
template<typename Type>
Matrix<complex<Type> > complexMatrix( const Matrix<Type>& );
template<typename Type>
Matrix<complex<Type> > complexMatrix( const Matrix<Type>&,
const Matrix<Type>& );
#include <matrix-impl.h>
}
// namespace splab
#endif
// MATRIX_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
*/
/*****************************************************************************
* matrix-impl.h
*
* Implementation for Matrix class.
*
* Zhang Ming, 2010-01 (revised 2010-12), Xi'an Jiaotong University.
*****************************************************************************/
/**
* initialize
*/
template <typename Type>
void Matrix<Type>::init( int rows, int columns )
{
nRow = rows;
nColumn = columns;
nTotal = nRow * nColumn;
pv0 = new Type[nTotal];
prow0 = new Type*[nRow];
prow1 = new Type*[nRow];
assert( pv0 != NULL );
assert( prow0 != NULL );
assert( prow1 != NULL );
Type *p = pv0;
pv1 = pv0 - 1;
for( int i=0; i<nRow; ++i )
{
prow0[i] = p;
prow1[i] = p-1;
p += nColumn;
}
prow1--;
}
/**
* copy matrix from normal array
*/
template <typename Type>
inline void Matrix<Type>::copyFromArray( const Type *v )
{
for( long i=0; i<nTotal; ++i )
pv0[i] = v[i];
}
/**
* set matrix by a scalar
*/
template <typename Type>
inline void Matrix<Type>::setByScalar( const Type &x )
{
for( long i=0; i<nTotal; ++i )
pv0[i] = x;
}
/**
* destroy the matrix
*/
template <typename Type>
void Matrix<Type>::destroy()
{
if( pv0 == NULL )
return ;
else
delete []pv0;
if( prow0 != NULL )
delete []prow0;
prow1++;
if( prow1 != NULL )
delete []prow1;
}
/**
* constructors and destructor
*/
template <typename Type>
Matrix<Type>::Matrix()
: pv0(0), pv1(0), prow0(0), prow1(0), nRow(0), nColumn(0), nTotal(0)
{
}
template <typename Type>
Matrix<Type>::Matrix( const Matrix<Type> &A )
{
init( A.nRow, A.nColumn );
copyFromArray( A.pv0 );
}
template <typename Type>
Matrix<Type>::Matrix( int rows, int columns, const Type &x )
{
init( rows,columns );
setByScalar(x);
}
template <typename Type>
Matrix<Type>::Matrix( int rows, int columns, const Type *arrays )
{
init( rows,columns );
copyFromArray( arrays );
}
template <typename Type>
Matrix<Type>::~Matrix()
{
destroy();
}
/**
* overload evaluate operator = from matrix to matrix
*/
template <typename Type>
Matrix<Type>& Matrix<Type>::operator=( const Matrix<Type> &A )
{
if( pv0 == A.pv0 )
return *this;
if( nRow == A.nRow && nColumn == A.nColumn )
copyFromArray( A.pv0 );
else
{
destroy();
init( A.nRow, A.nColumn );
copyFromArray( A.pv0 );
}
return *this;
}
/**
* overload evaluate operator = from scalar to matrix
*/
template <typename Type>
inline Matrix<Type>& Matrix<Type>::operator=( const Type &x )
{
setByScalar( x );
return *this;
}
/**
* overload operator [] for 0-offset access
*/
template <typename Type>
inline Type* Matrix<Type>::operator[]( int i )
{
#ifdef BOUNDS_CHECK
assert( 0 <= i );
assert( i < nRow );
#endif
return prow0[i];
}
template <typename Type>
inline const Type* Matrix<Type>::operator[]( int i ) const
{
#ifdef BOUNDS_CHECK
assert( 0 <= i );
assert( i < nRow );
#endif
return prow0[i];
}
/**
* overload operator () for 1-offset access
*/
template <typename Type>
inline Type& Matrix<Type>::operator()( int row, int column )
{
#ifdef BOUNDS_CHECK
assert( 1 <= row );
assert( row <= nRow ) ;
assert( 1 <= column);
assert( column <= nColumn );
#endif
return prow1[row][column];
}
template <typename Type>
inline const Type& Matrix<Type>::operator()( int row, int column ) const
{
#ifdef BOUNDS_CHECK
assert( 1 <= row );
assert( row <= nRow ) ;
assert( 1 <= column);
assert( column <= nColumn );
#endif
return prow1[row][column];
}
/**
* type conversion functions
*/
template <typename Type>
inline Matrix<Type>::operator Type*()
{
return pv0;
}
template <typename Type>
inline Matrix<Type>::operator const Type*() const
{
return pv0;
}
//template <typename Type>
//inline Matrix<Type>::operator Type**()
//{
// return prow0;
//}
//
//template <typename Type>
//inline Matrix<Type>::operator const Type**() const
//{
// return prow0;
//}
/**
* get the matrix's size
*/
template <typename Type>
inline long Matrix<Type>::size() const
{
return nTotal;
}
/**
* get the matrix's dimension
*/
template <typename Type>
int Matrix<Type>::dim( int dimension ) const
{
#ifdef BOUNDS_CHECK
assert( dimension >= 1);
assert( dimension <= 2);
#endif
if( dimension == 1 )
return nRow;
else if( dimension == 2 )
return nColumn;
else
return 0;
}
template <typename Type>
inline int Matrix<Type>::rows() const
{
return nRow;
}
template <typename Type>
inline int Matrix<Type>::cols() const
{
return nColumn;
}
/**
* reallocate matrix's size
*/
template <typename Type>
Matrix<Type>& Matrix<Type>::resize( int rows, int columns )
{
if( rows == nRow && columns == nColumn )
return *this;
destroy();
init( rows, columns );
return *this;
}
/**
* get the matrix's row vector
*/
template <typename Type>
Vector<Type> Matrix<Type>::getRow( int row ) const
{
#ifdef BOUNDS_CHECK
assert( row >= 0 );
assert( row < nRow );
#endif
Vector<Type> tmp( nColumn );
for( int j=0; j<nColumn; ++j )
tmp[j] = prow0[row][j];
return tmp;
}
/**
* get the matrix's column vector
*/
template <typename Type>
Vector<Type> Matrix<Type>::getColumn( int column ) const
{
#ifdef BOUNDS_CHECK
assert( column >= 0 );
assert( column < nColumn );
#endif
Vector<Type> tmp( nRow );
for( int i=0; i<nRow; ++i )
tmp[i] = prow0[i][column];
return tmp;
}
/**
* set the matrix's row vector
*/
template <typename Type>
void Matrix<Type>::setRow( const Vector<Type> &v, int row )
{
#ifdef BOUNDS_CHECK
assert( row >= 0 );
assert( row < nRow );
assert( v.dim() == nColumn );
#endif
for( int j=0; j<nColumn; ++j )
prow0[row][j] = v[j];
}
/**
* set the matrix's column vector
*/
template <typename Type>
void Matrix<Type>::setColumn( const Vector<Type> &v, int column )
{
#ifdef BOUNDS_CHECK
assert( column >= 0 );
assert( column < nColumn );
assert( v.dim() == nRow );
#endif
for( int i=0; i<nRow; ++i )
prow0[i][column] = v[i];
}
/**
* compound assignment operators +=
*/
template <typename Type>
Matrix<Type>& Matrix<Type>::operator+=( const Type &x )
{
Type **rowPtr = prow0;
Type *colPtr = 0;
for( int i=0; i<nRow; ++i )
{
colPtr = *rowPtr++;
for( int j=0; j<nColumn; ++j )
*colPtr++ += x;
}
return *this;
}
template <typename Type>
Matrix<Type>& Matrix<Type>::operator+=( const Matrix<Type> &rhs )
{
assert( nRow == rhs.rows() );
assert( nColumn == rhs.cols() );
Type **rowPtrL = prow0;
Type *colPtrL = 0;
Type **rowPtrR = rhs.prow0;
const Type *colPtrR = 0;
for( int i=0; i<nRow; ++i )
{
colPtrL = *rowPtrL++;
colPtrR = *rowPtrR++;
for( int j=0; j<nColumn; ++j )
*colPtrL++ += *colPtrR++;
}
return *this;
}
/**
* compound assignment operators -=
*/
template <typename Type>
Matrix<Type>& Matrix<Type>::operator-=( const Type &x )
{
Type **rowPtr = prow0;
Type *colPtr = 0;
for( int i=0; i<nRow; ++i )
{
colPtr = *rowPtr++;
for( int j=0; j<nColumn; ++j )
*colPtr++ -= x;
}
return *this;
}
template <typename Type>
Matrix<Type>& Matrix<Type>::operator-=( const Matrix<Type> &rhs )
{
assert( nRow == rhs.rows() );
assert( nColumn == rhs.cols() );
Type **rowPtrL = prow0;
Type *colPtrL = 0;
Type **rowPtrR = rhs.prow0;
const Type *colPtrR = 0;
for( int i=0; i<nRow; ++i )
{
colPtrL = *rowPtrL++;
colPtrR = *rowPtrR++;
for( int j=0; j<nColumn; ++j )
*colPtrL++ -= *colPtrR++;
}
return *this;
}
/**
* compound assignment operators *=
*/
template <typename Type>
Matrix<Type>& Matrix<Type>::operator*=( const Type &x )
{
Type **rowPtr = prow0;
Type *colPtr = 0;
for( int i=0; i<nRow; ++i )
{
colPtr = *rowPtr++;
for( int j=0; j<nColumn; ++j )
*colPtr++ *= x;
}
return *this;
}
// WARNING: this is element-by-element multiplication
template <typename Type>
Matrix<Type>& Matrix<Type>::operator*=( const Matrix<Type> &rhs )
{
assert( nRow == rhs.rows() );
assert( nColumn == rhs.cols() );
Type **rowPtrL = prow0;
Type *colPtrL = 0;
Type **rowPtrR = rhs.prow0;
const Type *colPtrR = 0;
for( int i=0; i<nRow; ++i )
{
colPtrL = *rowPtrL++;
colPtrR = *rowPtrR++;
for( int j=0; j<nColumn; ++j )
*colPtrL++ *= *colPtrR++;
}
return *this;
}
/**
* compound assignment operators /=
*/
template <typename Type>
Matrix<Type>& Matrix<Type>::operator/=( const Type &x )
{
Type **rowPtr = prow0;
Type *colPtr = 0;
for( int i=0; i<nRow; ++i )
{
colPtr = *rowPtr++;
for( int j=0; j<nColumn; ++j )
*colPtr++ /= x;
}
return *this;
}
// WARNING: this is element-by-element division
template <typename Type>
Matrix<Type>& Matrix<Type>::operator/=( const Matrix<Type> &rhs )
{
assert( nRow == rhs.rows() );
assert( nColumn == rhs.cols() );
Type **rowPtrL = prow0;
Type *colPtrL = 0;
Type **rowPtrR = rhs.prow0;
const Type *colPtrR = 0;
for( int i=0; i<nRow; ++i )
{
colPtrL = *rowPtrL++;
colPtrR = *rowPtrR++;
for( int j=0; j<nColumn; ++j )
*colPtrL++ /= *colPtrR++;
}
return *this;
}
/**
* Overload the output stream function.
*/
template <typename Type>
ostream& operator<<( ostream &out, const Matrix<Type> &A )
{
int rows = A.rows();
int columns = A.cols();
out << "size: " << rows << " by " << columns << "\n";
for( int i=0; i<rows; ++i )
{
for( int j=0; j<columns; ++j )
out << A[i][j] << "\t";
out << "\n";
}
return out;
}
/**
* Overload the intput stream function.
*/
template <typename Type>
istream& operator>>( istream &in, Matrix<Type> &A )
{
int rows, columns;
in >> rows >> columns;
if( !( rows == A.rows() && columns == A.cols() ) )
A.resize( rows, columns );
for( int i=0; i<rows; ++i )
for( int j=0; j<columns; ++j )
in >> A[i][j];
return in;
}
/**
* get negative matrix
*/
template<typename Type>
Matrix<Type> operator-( const Matrix<Type> &A )
{
int rows = A.rows();
int columns = A.cols();
Matrix<Type> tmp( rows, columns );
for( int i=0; i<rows; ++i )
for( int j=0; j<columns; ++j )
tmp[i][j] = -A[i][j];
return tmp;
}
/**
* matrix-scalar addition
*/
template<typename Type>
inline Matrix<Type> operator+( const Matrix<Type> &A, const Type &x )
{
Matrix<Type> tmp( A );
return tmp += x;
}
template<typename Type>
inline Matrix<Type> operator+( const Type &x, const Matrix<Type> &A )
{
return A + x;
}
/**
* matrix-matrix addition
*/
template<typename Type>
inline Matrix<Type> operator+( const Matrix<Type> &A1, const Matrix<Type> &A2 )
{
Matrix<Type> tmp( A1 );
return tmp += A2;
}
/**
* matrix-scalar subtraction
*/
template<typename Type>
inline Matrix<Type> operator-( const Matrix<Type> &A, const Type &x )
{
Matrix<Type> tmp( A );
return tmp -= x;
}
template<typename Type>
inline Matrix<Type> operator-( const Type &x, const Matrix<Type> &A )
{
Matrix<Type> tmp( A );
return -tmp += x;
}
/**
* matrix-matrix subtraction
*/
template<typename Type>
inline Matrix<Type> operator-( const Matrix<Type> &A1, const Matrix<Type> &A2 )
{
Matrix<Type> tmp( A1 );
return tmp -= A2;
}
/**
* matrix-scaling multiplication
*/
template <typename Type>
inline Matrix<Type> operator*( const Matrix<Type> &A, const Type &x )
{
Matrix<Type> tmp( A );
return tmp *= x;
}
template <typename Type>
inline Matrix<Type> operator*( const Type &x, const Matrix<Type> &A )
{
return A * x;
}
/**
* matrix-matrix multiplication
*/
template <typename Type>
Matrix<Type> operator*( const Matrix<Type> &A1, const Matrix<Type> &A2 )
{
assert( A1.cols() == A2.rows() );
int rows = A1.rows();
int columns = A2.cols();
// int K = A1.cols();
Matrix<Type> tmp( rows, columns );
// for( int i=0; i<rows; ++i )
// for( int j=0; j<columns; ++j )
// {
// tmp[i][j] = 0;
// for( int k=0; k<K; ++k )
// tmp[i][j] += A1[i][k] * A2[k][j];
// }
mult( A1, A2, tmp );
return tmp;
}
/**
* matrix-vector multiplication
*/
template <typename Type>
Vector<Type> operator*( const Matrix<Type> &A, const Vector<Type> &b )
{
assert( A.cols() == b.dim() );
int rows = A.rows();
// int columns = A.cols();
Vector<Type> tmp(rows);
// for( int i=0; i<rows; ++i )
// {
// Type sum = 0;
// for( int j=0; j<columns; ++j )
// sum += A[i][j] * v[j];
// tmp[i] = sum;
// }
mult( A, b, tmp );
return tmp;
}
/**
* matrix-scalar division
*/
template <typename Type>
inline Matrix<Type> operator/( const Matrix<Type> &A, const Type &x )
{
Matrix<Type> tmp( A );
return tmp /= x;
}
template <typename Type>
Matrix<Type> operator/( const Type &x, const Matrix<Type> &A )
{
int rows = A.rows();
int clumns = A.cols();
Matrix<Type> tmp( rows,clumns );
for( int i=0; i<rows; ++i )
for( int j=0; j<clumns; ++j )
tmp[i][j] = x / A[i][j];
return tmp;
}
/**
* This is an optimized version of matrix multiplication,
* where the destination matrix has already been allocated.
*/
template <typename Type>
Matrix<Type>& mult( const Matrix<Type> &A, const Matrix<Type> &B,
Matrix<Type> &C )
{
int M = A.rows();
int N = B.cols();
int K = A.cols();
assert( B.rows() == K );
C.resize( M, N );
Type sum;
const Type *pRow,
*pCol;
for( int i=0; i<M; i++ )
for( int j=0; j<N; ++j )
{
pRow = &A[i][0];
pCol = &B[0][j];
sum = 0;
for( int k=0; k<K; ++k )
{
sum += (*pRow) * (*pCol);
pRow++;
pCol += N;
}
C[i][j] = sum;
}
return C;
}
/**
* This is an optimized version of matrix and vector multiplication,
* where the destination vector has already been allocated.
*/
template <typename Type>
Vector<Type>& mult( const Matrix<Type> &A, const Vector<Type> &b,
Vector<Type> &c )
{
int M = A.rows();
int N = A.cols();
assert( b.size() == N );
c.resize( M );
Type sum;
const Type *pRow,
*pCol;
for( int i=0; i<M; i++ )
{
pRow = &A[i][0];
pCol = &b[0];
sum = 0;
for( int j=0; j<N; ++j )
{
sum += (*pRow) * (*pCol);
pRow++;
pCol++;
}
c[i] = sum;
}
return c;
}
/**
* matrix-matrix elementwise multiplication
*/
template<typename Type>
inline Matrix<Type> elemMult( const Matrix<Type> &A1, const Matrix<Type> &A2 )
{
Matrix<Type> tmp( A1 );
return tmp *= A2;
}
template <typename Type>
inline Matrix<Type>& elemMultEq( Matrix<Type> &A1, const Matrix<Type> &A2 )
{
return A1 *= A2;
}
/**
* matrix-matrix elementwise division
*/
template <typename Type>
inline Matrix<Type> elemDivd( const Matrix<Type> &A1, const Matrix<Type> &A2 )
{
Matrix<Type> tmp( A1 );
return tmp /= A2;
}
template <typename Type>
inline Matrix<Type>& elemDivdEq( Matrix<Type> &A1, const Matrix<Type> &A2 )
{
return A1 /= A2;
}
/**
* matrix tranpose
*/
template <typename Type>
Matrix<Type> trT( const Matrix<Type> &A )
{
int rows = A.cols();
int clumns = A.rows();
Matrix<Type> tmp( rows, clumns );
for( int i=0; i<rows; ++i )
for( int j=0; j<clumns; ++j )
tmp[i][j] = A[j][i];
return tmp;
}
/**
* matrix conjugate tranpose
*/
template <typename Type>
Matrix<Type> trH( const Matrix<Type> &A )
{
int rows = A.cols();
int clumns = A.rows();
Matrix<Type> tmp( rows, clumns );
for( int i=0; i<rows; ++i )
for( int j=0; j<clumns; ++j )
tmp[i][j] = conj(A[j][i]);
return tmp;
}
/**
* matrix-matrix tranpose multiplication: A^T * B.
*/
template <typename Type>
Matrix<Type> trMult( const Matrix<Type> &A1, const Matrix<Type> &A2 )
{
assert( A1.rows() == A2.rows() );
int rows = A1.cols();
int columns = A2.cols();
int K = A1.rows();
Matrix<Type> tmp( rows, columns );
// for( int i=0; i<rows; ++i )
// for( int j=0; j<columns; ++j )
// {
// Type sum = 0;
// for( int k=0; k<K; ++k )
// sum += A1[k][i] * A2[k][j];
// tmp[i][j] = sum;
// }
for( int i=0; i<rows; ++i )
for( int j=0; j<columns; ++j )
for( int k=0; k<K; ++k )
tmp[i][j] += A1[k][i] * A2[k][j];
return tmp;
}
/**
* matrix-vector tranpose multiplication: A^T * b.
*/
template <typename Type>
Vector<Type> trMult( const Matrix<Type> &A, const Vector<Type> &v )
{
assert( A.rows() == v.dim() );
int rows = A.rows();
int columns = A.cols();
Vector<Type> tmp( columns );
// for( int i=0; i<columns; ++i )
// {
// Type sum = 0;
// for( int j=0; j<rows; ++j )
// sum += A[j][i] * v[j];
// tmp[i] = sum;
// }
for( int i=0; i<columns; ++i )
for( int j=0; j<rows; ++j )
tmp[i] += A[j][i] * v[j];
return tmp;
}
/**
* matrix-matrix tranpose multiplication: A * B^T.
*/
template <typename Type>
Matrix<Type> multTr( const Matrix<Type> &A1, const Matrix<Type> &A2 )
{
assert( A1.cols() == A2.cols() );
int rows = A1.rows();
int columns = A2.rows();
int K = A1.cols();
Matrix<Type> tmp( rows, columns );
// for( int i=0; i<rows; ++i )
// for( int j=0; j<columns; ++j )
// {
// Type sum = 0;
// for( int k=0; k<K; ++k )
// sum += A1[i][k] * A2[j][k];
// tmp[i][j] = sum;
// }
for( int i=0; i<rows; ++i )
for( int j=0; j<columns; ++j )
for( int k=0; k<K; ++k )
tmp[i][j] += A1[i][k] * A2[j][k];
return tmp;
}
/**
* vector-vector tranpose multiplication: a * b^T.
*/
template <typename Type>
Matrix<Type> multTr( const Vector<Type> &a, const Vector<Type> &b )
{
int rows = a.dim();
int columns = b.dim();
Matrix<Type> tmp( rows, columns );
for( int i=0; i<rows; ++i )
for( int j=0; j<columns; ++j )
tmp[i][j] = a[i]*b[j];
return tmp;
}
/**
* matrix-matrix tranpose multiplication: A^H * B.
*/
template <typename Type>
Matrix<complex<Type> > trMult( const Matrix<complex<Type> > &A1,
const Matrix<complex<Type> > &A2 )
{
assert( A1.rows() == A2.rows() );
int rows = A1.cols();
int columns = A2.cols();
int K = A1.rows();
Matrix<complex<Type> > tmp( rows, columns );
// for( int i=0; i<rows; ++i )
// for( int j=0; j<columns; ++j )
// {
// Type sum = 0;
// for( int k=0; k<K; ++k )
// sum += A1[k][i] * A2[k][j];
// tmp[i][j] = sum;
// }
for( int i=0; i<rows; ++i )
for( int j=0; j<columns; ++j )
for( int k=0; k<K; ++k )
tmp[i][j] += conj(A1[k][i]) * A2[k][j];
return tmp;
}
/**
* matrix-vector tranpose multiplication: A^H * b.
*/
template <typename Type>
Vector<complex<Type> > trMult( const Matrix<complex<Type> > &A,
const Vector<complex<Type> > &v )
{
assert( A.rows() == v.dim() );
int rows = A.rows();
int columns = A.cols();
Vector<complex<Type> > tmp( columns );
// for( int i=0; i<columns; ++i )
// {
// Type sum = 0;
// for( int j=0; j<rows; ++j )
// sum += A[j][i] * v[j];
// tmp[i] = sum;
// }
for( int i=0; i<columns; ++i )
for( int j=0; j<rows; ++j )
tmp[i] += conj(A[j][i]) * v[j];
return tmp;
}
/**
* matrix-matrix tranpose multiplication: A * B^H.
*/
template <typename Type>
Matrix<complex<Type> > multTr( const Matrix<complex<Type> > &A1,
const Matrix<complex<Type> > &A2 )
{
assert( A1.cols() == A2.cols() );
int rows = A1.rows();
int columns = A2.rows();
int K = A1.cols();
Matrix<complex<Type> > tmp( rows, columns );
// for( int i=0; i<rows; ++i )
// for( int j=0; j<columns; ++j )
// {
// Type sum = 0;
// for( int k=0; k<K; ++k )
// sum += A1[i][k] * A2[j][k];
// tmp[i][j] = sum;
// }
for( int i=0; i<rows; ++i )
for( int j=0; j<columns; ++j )
for( int k=0; k<K; ++k )
tmp[i][j] += A1[i][k] * conj(A2[j][k]);
return tmp;
}
/**
* vector-vector tranpose multiplication: a * b^H.
*/
template <typename Type>
Matrix<complex<Type> > multTr( const Vector<complex<Type> > &a,
const Vector<complex<Type> > &b )
{
int rows = a.dim();
int columns = b.dim();
Matrix<complex<Type> > tmp( rows, columns );
for( int i=0; i<rows; ++i )
for( int j=0; j<columns; ++j )
tmp[i][j] = a[i]*conj(b[j]);
return tmp;
}
/**
* Generate the identity matrix.
*/
template <typename Type>
Matrix<Type> eye( int N, const Type &x )
{
Matrix<Type> tmp( N, N );
for( int i=0; i<N; ++i )
tmp[i][i] = x;
return tmp;
}
/**
* Get the diagonal entries of matrix.
*/
template <typename Type>
Vector<Type> diag( const Matrix<Type> &A )
{
int nColumn = A.rows();
if( nColumn > A.cols() )
nColumn = A.cols();
Vector<Type> tmp( nColumn );
for( int i=0; i<nColumn; ++i )
tmp[i] = A[i][i];
return tmp;
}
/**
* Generate the diagonal of matrix by given its diagonal elements.
*/
template <typename Type>
Matrix<Type> diag( const Vector<Type> &d )
{
int N = d.size();
Matrix<Type> tmp( N, N );
for( int i=0; i<N; ++i )
tmp[i][i] = d[i];
return tmp;
}
/**
* Compute Frobenius norm of matrix.
*/
template <typename Type>
Type norm( const Matrix<Type> &A )
{
int m = A.rows();
int n = A.cols();
Type sum = 0;
for( int i=1; i<=m; ++i )
for( int j=1; j<=n; ++j )
sum += A(i,j) * A(i,j);
return sqrt(sum);
}
template <typename Type>
Type norm( const Matrix<complex<Type> > &A )
{
int m = A.rows();
int n = A.cols();
Type sum = 0;
for( int i=1; i<=m; ++i )
for( int j=1; j<=n; ++j )
sum += norm(A(i,j));
return sqrt(sum);
}
/**
* Swap two matrixes.
*/
template <typename Type> void swap( Matrix<Type> &lhs, Matrix<Type> &rhs )
{
int m = lhs.rows();
int n = lhs.cols();
assert( m == rhs.rows() );
assert( n == rhs.cols() );
for( int i=1; i<=m; ++i )
for( int j=1; j<=n; ++j )
swap( lhs(i,j), rhs(i,j) );
}
/**
* Matrix's column vecotrs sum.
*/
template <typename Type>
Vector<Type> sum( const Matrix<Type> &A )
{
int m = A.rows();
int n = A.cols();
Vector<Type> sum(n);
for( int j=1; j<=n; ++j )
for( int i=1; i<=m; ++i )
sum(j) += A(i,j);
return sum;
}
/**
* Minimum of matrix's column vecotrs.
*/
template <typename Type>
Vector<Type> min( const Matrix<Type> &A )
{
int m = A.rows();
int n = A.cols();
Vector<Type> sum(n);
for( int j=1; j<=n; ++j )
{
Type tmp = A(1,j);
for( int i=2; i<m; ++i )
if( tmp > A(i,j) )
tmp = A(i,j);
sum(j) = tmp;
}
return sum;
}
/**
* Maximum of matrix's column vecotrs.
*/
template <typename Type>
Vector<Type> max( const Matrix<Type> &A )
{
int m = A.rows();
int n = A.cols();
Vector<Type> sum(n);
for( int j=1; j<=n; ++j )
{
Type tmp = A(1,j);
for( int i=2; i<m; ++i )
if( tmp < A(i,j) )
tmp = A(i,j);
sum(j) = tmp;
}
return sum;
}
/**
* Matrix's column vecotrs mean.
*/
template <typename Type>
inline Vector<Type> mean( const Matrix<Type> &A )
{
return sum(A) / Type(A.rows());
}
/**
* Convert real matrix to complex matrix.
*/
template <typename Type>
Matrix<complex<Type> > complexMatrix( const Matrix<Type> &rA )
{
int rows = rA.rows();
int columns = rA.cols();
Matrix<complex<Type> > cA( rows, columns );
for( int i=0; i<rows; ++i )
for( int j=0; j<columns; ++j )
cA[i][j] = rA[i][j];
return cA;
}
template <typename Type>
Matrix<complex<Type> > complexMatrix( const Matrix<Type> &mR,
const Matrix<Type> &mI )
{
int rows = mR.rows();
int columns = mR.cols();
assert( rows == mI.rows() );
assert( columns == mI.cols() );
Matrix<complex<Type> > cA( rows, columns );
for( int i=0; i<rows; ++i )
for( int j=0; j<columns; ++j )
cA[i][j] = complex<Type>( mR[i][j], mI[i][j] );
return cA;
}
/**
* Get magnitude of a complex matrix.
*/
template <typename Type>
Matrix<Type> abs( const Matrix<complex<Type> > &A )
{
int m = A.rows(),
n = A.cols();
Matrix<Type> tmp( m, n );
for( int i=0; i<m; ++i )
for( int j=0; j<n; ++j )
tmp[i][j] = abs( A[i][j] );
return tmp;
}
/**
* Get angle of a complex matrix.
*/
template <typename Type>
Matrix<Type> arg( const Matrix<complex<Type> > &A )
{
int m = A.rows(),
n = A.cols();
Matrix<Type> tmp( m, n );
for( int i=0; i<m; ++i )
for( int j=0; j<n; ++j )
tmp[i][j] = arg( A[i][j] );
return tmp;
}
/**
* Get real part of a complex matrix.
*/
template <typename Type>
Matrix<Type> real( const Matrix<complex<Type> > &A )
{
int m = A.rows(),
n = A.cols();
Matrix<Type> tmp( m, n );
for( int i=0; i<m; ++i )
for( int j=0; j<n; ++j )
tmp[i][j] = A[i][j].real();
return tmp;
}
/**
* Get imaginary part of a complex matrix.
*/
template <typename Type>
Matrix<Type> imag( const Matrix<complex<Type> > &A )
{
int m = A.rows(),
n = A.cols();
Matrix<Type> tmp( m, n );
for( int i=0; i<m; ++i )
for( int j=0; j<n; ++j )
tmp[i][j] = A[i][j].imag();
return tmp;
}
测试代码:
/*****************************************************************************
* matrix_test.cpp
*
* Matrix class testing.
*
* Zhang Ming, 2010-01 (revised 2010-12), Xi'an Jiaotong University.
*****************************************************************************/
#define BOUNDS_CHECK
#include <iostream>
#include <iomanip>
#include <matrix.h>
using namespace std;
using namespace splab;
typedef double Type;
const int M = 3;
const int N = 3;
int main()
{
Type x;
Matrix<Type> m1;
m1.resize( M, N );
x = 1.0;
m1 = x;
cout << "matrix m1 : " << m1 << endl;
x = 2.0;
Matrix<Type> m2( M, N, x );
cout << "matrix m2 : " << m2 << endl;
Matrix<Type> m3 = m1;
cout << "matrix m3 : " << m3 << endl;
m3.resize( 3, 4 );
for( int i=1; i<=3; ++i )
for( int j=1; j<=4; ++j )
m3(i,j) = Type(i*j);
int row = m3.dim(1);
int column = m3.dim(2);
cout << "the row number of new matrix m3 : " << row << endl;
cout << "the column number of new matrix m3 : " << column << endl;
cout << "new matrix m3 : " << m3 << endl;
cout << "the diagonal matrix of m3 : " << diag( m3 ) << endl;
cout << "the transpose matrix of m3 : " << trT( m3 ) << endl;
cout << endl << "\t\t\t\tmatrix-scalar operand" << endl << endl;
cout << "scalar x = " << x << endl;
cout << "m1 + x : " << m1+x << endl;
cout << "x + m1 : " << x+m1 << endl;
m1 += x;
cout << "m1 += x : " << m1 << endl;
cout << "m1 - x : " << m1-x << endl;
cout << "x - m1 : " << x-m1 << endl;
m1 -= x;
cout << "m1 -= x : " << m1 << endl;
cout << "m1 * x : " << m1*x << endl;
cout << "x * m1 : " << x*m1 << endl;
m1 *= x;
cout << "m1 *= x : " << m1 << endl;
cout << "m1 / x : " << m1/x << endl;
cout << "x / m1 : " << x/m1 << endl;
m1 /= x;
cout << "m1 /= x : " << m1 << endl;
cout << endl << "\t\t\telementwise matrix-matrix operand" << endl << endl;
cout << "m1 + m2 : " << m1 + m2 << endl;
m1 += m2;
cout << "m1 += m2 : " << m1 << endl;
cout << "m1 - m2 : " << m1-m2 << endl;
m1 -= m2;
cout << "m1 -= m2 : " << m1 << endl;
m1 *= m2;
cout << "m1 *= m2 : " << m1 << endl;
m1 /= m2;
cout << "m1 /= m2 : " << m1 << endl;
cout << "column minimum vector of m1 : " << min(m1) << endl;
cout << "column maximum vector of m1 : " << max(m1) << endl;
cout << "column sum vector of m1 : " << sum(m1) << endl;
cout << "column mean vector of m1 : " << mean(m1) << endl;
cout << endl << "\t\t\t\tmatrix-vector operand" << endl << endl;
Vector<Type> v1( 3, 2 );
cout << "vector v1 : " << v1 << endl;
cout << "m1 * v1 : " << m1*v1 << endl;
cout << "m1^T * v1 : " << trMult(m1, v1) << endl;
cout << endl << "\t\t\t\tmatrix-matrix operand" << endl << endl;
cout << "m1 * m2 : " << m1*m2 << endl;
cout << "m3^T * m2 : " << trMult(m3, m2) << endl;
cout << "m1 * m3^T^T : " << multTr(m1, trT(m3)) << endl;
Matrix<int> m4( 4, 5 );
Vector<int> v2(5);
v2[0] = 1; v2[1] = 2; v2[2] = 3; v2[3] = 4; v2[4] = 5;
for( int i=0; i<4; ++i )
m4.setRow( i*v2, i );
cout << "matrix m4 : " << m4 << endl;
cout << "column vectors of m4 : " << endl;
for( int j=0; j<5; ++j )
cout << "the " << j << "th column" << m4.getColumn(j) << endl;
v2.resize(4);
v2[0] = 1; v2[1] = 2; v2[2] = 3; v2[3] = 4;
for( int j=0; j<5; ++j )
m4.setColumn( j*v2, j );
cout << "row vectors of m4 : " << endl;
for( int i=0; i<4; ++i )
cout << "the " << i << "th row" << m4.getRow(i) << endl;
cout << endl << "\t\t\t\tcomplex matrix operand" << endl << endl;
complex<Type> c = polar(1.0,PI/4);
Matrix< complex<Type> > A( M, N ), B(M,N), C( M, N+1 );
for( int i=0; i<M; i++ )
for( int j=0; j<N; j++ )
{
A[i][j] = complex<Type>( Type(0.3*i+0.7*j), sin(Type(i+j)) );
B[i][j] = complex<Type>( cos(Type(i+j)), Type(0.8*i+0.2*j) );
}
for( int i=0; i<N; ++i )
C.setColumn( A.getColumn(i), i );
C.setColumn( B.getColumn(0), N );
cout << setiosflags(ios::fixed) << setprecision(2);
cout << "Matrix A: " << A << endl;
cout << "Matrix B: " << B << endl;
cout << "Matrix C: " << C << endl;
cout << "Absolute of C : " << abs(C) << endl;
cout << "Angle of C : " << arg(C) << endl;
cout << "Real part of C : " << real(C) << endl;
cout << "Imaginary part of C : " << imag(C) << endl;
cout << "c*A - A*B + trH(B) " << c*A - A*B + trH(B) << endl;
cout << "A*C * C^H*B " << A*C * trMult(C,B) << endl;
cout << "(A.*B+C*C^H)) ./ A^T " << elemDivd( elemMult(A,B)+multTr(C,C), trT(A) ) << endl;
cout << "diag( A * diag(C) ) " << diag( A*diag(C)/c ) * c << endl;
cout << "A .* ( B./(c*m1) ) " << elemMultEq( A, elemDivdEq(B,c*complexMatrix(m1)) ) << endl;
return 0;
}
运行结果:
3 6 9 12
the diagonal matrix of m3 : size: 3 by 1
1
4
9
the transpose matrix of m3 : size: 4 by 3
1 2 3
2 4 6
3 6 9
4 8 12
matrix-scalar operand
scalar x = 2
m1 + x : size: 3 by 3
3 3 3
3 3 3
3 3 3
x + m1 : size: 3 by 3
3 3 3
3 3 3
3 3 3
m1 += x : size: 3 by 3
3 3 3
3 3 3
3 3 3
m1 - x : size: 3 by 3
1 1 1
1 1 1
1 1 1
x - m1 : size: 3 by 3
-1 -1 -1
-1 -1 -1
-1 -1 -1
m1 -= x : size: 3 by 3
1 1 1
1 1 1
1 1 1
m1 * x : size: 3 by 3
2 2 2
2 2 2
2 2 2
x * m1 : size: 3 by 3
2 2 2
2 2 2
2 2 2
m1 *= x : size: 3 by 3
2 2 2
2 2 2
2 2 2
m1 / x : size: 3 by 3
1 1 1
1 1 1
1 1 1
x / m1 : size: 3 by 3
1 1 1
1 1 1
1 1 1
m1 /= x : size: 3 by 3
1 1 1
1 1 1
1 1 1
elementwise matrix-matrix operand
m1 + m2 : size: 3 by 3
3 3 3
3 3 3
3 3 3
m1 += m2 : size: 3 by 3
3 3 3
3 3 3
3 3 3
m1 - m2 : size: 3 by 3
1 1 1
1 1 1
1 1 1
m1 -= m2 : size: 3 by 3
1 1 1
1 1 1
1 1 1
m1 *= m2 : size: 3 by 3
2 2 2
2 2 2
2 2 2
m1 /= m2 : size: 3 by 3
1 1 1
1 1 1
1 1 1
column minimum vector of m1 : size: 3 by 1
1
1
1
column maximum vector of m1 : size: 3 by 1
1
1
1
column sum vector of m1 : size: 3 by 1
3
3
3
column mean vector of m1 : size: 3 by 1
1
1
1
matrix-vector operand
vector v1 : size: 3 by 1
2
2
2
m1 * v1 : size: 3 by 1
6
6
6
m1^T * v1 : size: 3 by 1
6
6
6
matrix-matrix operand
m1 * m2 : size: 3 by 3
6 6 6
6 6 6
6 6 6
m3^T * m2 : size: 4 by 3
12 12 12
24 24 24
36 36 36
48 48 48
m1 * m3^T^T : size: 3 by 4
6 12 18 24
6 12 18 24
6 12 18 24
matrix m4 : size: 4 by 5
0 0 0 0 0
1 2 3 4 5
2 4 6 8 10
3 6 9 12 15
column vectors of m4 :
the 0th columnsize: 4 by 1
0
1
2
3
the 1th columnsize: 4 by 1
0
2
4
6
the 2th columnsize: 4 by 1
0
3
6
9
the 3th columnsize: 4 by 1
0
4
8
12
the 4th columnsize: 4 by 1
0
5
10
15
row vectors of m4 :
the 0th rowsize: 5 by 1
0
1
2
3
4
the 1th rowsize: 5 by 1
0
2
4
6
8
the 2th rowsize: 5 by 1
0
3
6
9
12
the 3th rowsize: 5 by 1
0
4
8
12
16
complex matrix operand
Matrix A: size: 3 by 3
(0.00,0.00) (0.70,0.84) (1.40,0.91)
(0.30,0.84) (1.00,0.91) (1.70,0.14)
(0.60,0.91) (1.30,0.14) (2.00,-0.76)
Matrix B: size: 3 by 3
(1.00,0.00) (0.54,0.20) (-0.42,0.40)
(0.54,0.80) (-0.42,1.00) (-0.99,1.20)
(-0.42,1.60) (-0.99,1.80) (-0.65,2.00)
Matrix C: size: 3 by 4
(0.00,0.00) (0.70,0.84) (1.40,0.91) (1.00,0.00)
(0.30,0.84) (1.00,0.91) (1.70,0.14) (0.54,0.80)
(0.60,0.91) (1.30,0.14) (2.00,-0.76) (-0.42,1.60)
Absolute of C : size: 3 by 4
0.00 1.09 1.67 1.00
0.89 1.35 1.71 0.97
1.09 1.31 2.14 1.65
Angle of C : size: 3 by 4
0.00 0.88 0.58 0.00
1.23 0.74 0.08 0.98
0.99 0.11 -0.36 1.83
Real part of C : size: 3 by 4
0.00 0.70 1.40 1.00
0.30 1.00 1.70 0.54
0.60 1.30 2.00 -0.42
Imaginary part of C : size: 3 by 4
0.00 0.84 0.91 0.00
0.84 0.91 0.14 0.80
0.91 0.14 -0.76 1.60
c*A - A*B + trH(B) size: 3 by 3
(3.33,-2.88) (4.60,-1.68) (4.37,-2.18)
(0.98,-4.19) (2.92,-3.71) (4.05,-3.88)
(-2.20,-4.87) (0.99,-6.38) (3.16,-6.90)
A*C * C^H*B size: 3 by 3
(-22.45,41.39) (-57.34,30.67) (-75.11,36.86)
(-13.03,52.79) (-53.59,46.63) (-73.18,52.60)
(16.66,60.01) (-25.93,77.12) (-40.26,93.59)
(A.*B+C*C^H)) ./ A^T size: 3 by 3
(inf,nan) (3.20,-4.47) (2.78,-0.12)
(2.30,-2.85) (3.12,-2.22) (3.99,2.29)
(-0.04,-1.18) (3.05,0.07) (3.81,3.69)
diag( A * diag(C) ) size: 3 by 3
(3.42,2.24) (0.00,0.00) (0.00,0.00)
(0.00,0.00) (3.68,0.81) (0.00,0.00)
(0.00,0.00) (0.00,0.00) (4.60,-1.70)
A .* ( B./(c*m1) ) size: 3 by 3
(0.00,0.00) (0.57,0.27) (-0.54,0.80)
(0.13,0.85) (-0.50,1.38) (0.03,2.65)
(-0.79,1.62) (0.47,2.65) (3.32,3.03)
Process returned 0 (0x0) execution time : 0.281 s
Press any key to continue.