/home5/piconepr/public_html/isip/projects/speech/software/documentation/class/math/matrix/MMatrix/mmat_09.cc
// file: $isip/class/math/matrix/MMatrix/mmat_09.cc
// version: $Id: mmat_09.cc 7924 2002-03-15 23:30:38Z zheng $
//
// isip include files
//
#include "MMatrixMethods.h"
#include "MMatrix.h"
[snipped]
// method: sub
//
// arguments:
// MMatrix<TScalar, TIntegral>& this: (output) class operand
// const MMatrix<TScalar, TIntegral>& m1: (input) operand 1
// const MMatrix<TAScalar, TAIntegral>& m2: (input) operand 2
//
// return: a boolean value indicating status
//
// subtract the two argument matrices and store result in the current matrix
// for example:
//
// [1 2 3] [2 0 4] [-1 2 -1]
// m1 = [4 5 6] m2 = [0 6 0] m_out = m1 - m2 = [ 4 -1 6]
// [7 8 9] [8 0 1] [-1 8 8]
//
// this method has evolved somewhat as our ability to efficiently
// copy and assign matrices has improved. this version implements
// some complicated logic to do things as efficiently as possible.
// the approach can best be summarized by the following flowchart and matrix:
//
// (1) special cases:
// types are equal: process internal data directly for all types
// one of the types is sparse: process internal data directly
// one of the types is full: the result is always a full matrix
//
// (2) other cases:
// this matrix describes the type of the output matrix as a function
// of the input types. For the operation c = a - b,
// type1 is the type of a; type 2 is the type of b.
//
// TYPE2 --------> FULL DIAG SYMM LTRI UTRI SPRS
// TYPE1 -> FULL: full full full full full full
// DIAG: full diag symm ltri utri sprs
// SYMM: full symm symm full full full
// LTRI: full ltri full ltri full full
// UTRI: full utri full full utri full
// SPRS: full sprs full full full sprs
//
// the code that implements each combination attempts to exploit the
// properties of the matrices as much as possible.
//
template<class TScalar, class TIntegral, class TAScalar, class TAIntegral>
boolean MMatrixMethods::sub(MMatrix<TScalar, TIntegral>& this_a,
const MMatrix<TScalar, TIntegral>& m1_a,
const MMatrix<TAScalar, TAIntegral>& m2_a) {
// check if the two input matrices have the dimensions
//
if (!m1_a.checkDimensions(m2_a)) {
m1_a.debug(L"m1_a");
m2_a.debug(L"m2_a");
return Error::handle(name(), L"sub", Error::ARG, __FILE__, __LINE__);
}
// declare local variables
//
long nrows = m1_a.nrows_d;
long ncols = m2_a.ncols_d;
Integral::MTYPE old_type = this_a.type_d;
Integral::MTYPE type1 = m1_a.type_d;
Integral::MTYPE type2 = m2_a.type_d;
//---------------------------------------------------------------------------
// special cases:
// in this section, we deal with special cases that share the same
// code structure.
//---------------------------------------------------------------------------
// if the two types are the same, we can simply operate on the storage
// vectors directly
//
if (type1 == type2) {
// type: FULL, DIAGONAL, SYMMETRIC, LOWER_TRIANGULAR, UPPER_TRIANGULAR
//
if ((type1 == Integral::FULL) ||
(type1 == Integral::DIAGONAL) ||
(type1 == Integral::SYMMETRIC) ||
(type1 == Integral::LOWER_TRIANGULAR) ||
(type1 == Integral::UPPER_TRIANGULAR)) {
this_a.setDimensions(nrows, ncols, false, m1_a.type_d);
this_a.m_d.sub(m1_a.m_d, m2_a.m_d);
}
// type: SPARSE
//
else {
// set the output matrix to SPARSE
//
this_a.assign(m1_a, Integral::SPARSE);
// loop through each element of the second argument
//
for (long index = 0; index < m2_a.m_d.length(); index++) {
long row = m2_a.row_index_d(index);
long col = m2_a.col_index_d(index);
this_a.setValue(row, col, (TScalar)((TIntegral)(this_a.getValue(row, col) - m2_a.m_d(index))));
}
}
}
// if the first type is sparse, and the second type isn't sparse,
// we assign the output type to FULL (except for DIAG), copy the elements
// of the first, and then subtract elements of the second
//
else if ((type1 == Integral::SPARSE) && (type2 != Integral::SPARSE)) {
// create a FULL matrix for all types except DIAGONAL, and zero it
//
if (type2 == Integral::DIAGONAL) {
this_a.setDimensions(nrows, ncols, false, Integral::SPARSE);
}
else {
this_a.setDimensions(nrows, ncols, false, Integral::FULL);
}
this_a.clear(Integral::RETAIN);
// assign elements of the first sparse matrix
//
for (long index = 0; index < m1_a.m_d.length(); index++) {
long row = m1_a.row_index_d(index);
long col = m1_a.col_index_d(index);
this_a.setValue(row, col, m1_a.m_d(index));
}
// subtract elements of the second matrix
//
for (long row = 0; row < nrows; row++) {
// subtract the values
//
long start_col = m2_a.startColumn(row, this_a.type_d, type2);
long stop_col = m2_a.stopColumn(row, this_a.type_d, type2);
for (long col = start_col; col <= stop_col; col++) {
this_a.setValue(row, col, (TScalar)((TIntegral)(this_a(row, col) - m2_a(row, col))));
}
}
}
// if the second type is sparse, and the first type isn't sparse,
// we assign the output type to FULL (except for DIAGONAL), copy the
// elements of the first, and then subtract elements of the second
//
else if ((type1 != Integral::SPARSE) && (type2 == Integral::SPARSE)) {
// create a FULL matrix for all types except DIAGONAL, and zero it
//
if (m1_a.type_d == Integral::DIAGONAL) {
this_a.assign(m1_a, Integral::SPARSE);
}
else {
this_a.assign(m1_a, Integral::FULL);
}
// subtract elements of the second sparse matrix
//
for (long index = 0; index < m2_a.m_d.length(); index++) {
long row = m2_a.row_index_d(index);
long col = m2_a.col_index_d(index);
this_a.setValue(row, col, (TScalar)((TIntegral)(this_a(row, col) - m2_a(row, col))));
}
}
// if either type is FULL, we assign the output type to FULL,
// copy the elements from the first, then subtract elements of the second
//
else if (type1 == Integral::FULL) {
// create a FULL matrix through assignment
//
this_a.assign(m1_a, Integral::FULL);
// subtract elements of the second matrix
//
for (long row = 0; row < nrows; row++) {
// subtract the values
//
long start_col = m2_a.startColumn(row, this_a.type_d, type2);
long stop_col = m2_a.stopColumn(row, this_a.type_d, type2);
for (long col = start_col; col <= stop_col; col++) {
this_a.setValue(row, col, (TScalar)(TIntegral)((this_a(row, col) - m2_a(row, col))));
}
}
}
else if (type2 == Integral::FULL) {
// create a FULL matrix through assignment
//
this_a.assign(m1_a, Integral::FULL);
// subtract elements of the second matrix
//
for (long row = 0; row < nrows; row++) {
for (long col = 0; col < ncols; col++) {
this_a.setValue(row, col, (TScalar)((TIntegral)(this_a(row, col) - m2_a(row, col))));
}
}
}
//---------------------------------------------------------------------------
// in this section of the code, we deal with the remaining cases where
// the matrices are of different types, and some intelligent conversion
// must be done.
//---------------------------------------------------------------------------
// type1: DIAGONAL
// the remaining types are SYMMETRIC and *_TRIANGULAR.
// the output types follow the type of the second operand.
// the approach is to create a matrix, copy the diagonal values,
// and subtract the second operand.
//
else if (type1 == Integral::DIAGONAL) {
// set up the output matrix and assign the diagonal elements
//
this_a.setDimensions(nrows, ncols, false, type2);
this_a.clear(Integral::RETAIN);
this_a.setDiagonal(m1_a);
// (1) subtract the diagonal for all
//
for (long row = 0; row < nrows; row++) {
this_a.m_d(this_a.index(row, row)) -=
(TIntegral)m2_a.m_d(m2_a.index(row, row));
}
// (2) subtract the lower part for SYMMETRIC and LOWER
//
if ((type2 == Integral::SYMMETRIC) ||
(type2 == Integral::LOWER_TRIANGULAR)) {
for (long row = 0; row < nrows; row++) {
for (long col = 0; col < row; col++) {
this_a.m_d(this_a.index(row, col)) -=
(TIntegral)m2_a.m_d(m2_a.index(row, col));
}
}
}
// (3) subtract the upper part for SYMMETRIC and UPPER
//
if (type2 == Integral::UPPER_TRIANGULAR) {
for (long row = 0; row < nrows; row++) {
for (long col = row + 1; col < ncols; col++) {
this_a.m_d(this_a.index(row, col)) -=
(TIntegral)m2_a.m_d(m2_a.index(row, col));
}
}
}
}
// type2: SYMMETRIC, *_TRIANGULAR
// the remaining types are DIAGONAL, SYMMETRIC, and *_TRIANGULAR.
// the output types vary with the input.
// the approach is to copy the values of the first operand,
// and subtract the values of the second operand.
//
else {
// set up the output matrix: a diagonal matrix preserves the type
// of the input; all other types force a full matrix.
//
if ((type1 == Integral::SYMMETRIC) &&
(type2 == Integral::DIAGONAL)) {
this_a.assign(m1_a, Integral::SYMMETRIC);
}
if ((type1 == Integral::LOWER_TRIANGULAR) &&
(type2 == Integral::DIAGONAL)) {
this_a.assign(m1_a, Integral::LOWER_TRIANGULAR);
}
else if ((type1 == Integral::UPPER_TRIANGULAR) &&
(type2 == Integral::DIAGONAL)) {
this_a.assign(m1_a, Integral::UPPER_TRIANGULAR);
}
else {
this_a.assign(m1_a, Integral::FULL);
}
// subtract elements of the second matrix
//
for (long row = 0; row < nrows; row++) {
// subtract the values
//
long start_col = m2_a.startColumn(row, this_a.type_d, type2);
long stop_col = m2_a.stopColumn(row, this_a.type_d, type2);
for (long col = start_col; col <= stop_col; col++) {
this_a.m_d(this_a.index(row, col)) -=
(TIntegral)m2_a.m_d(m2_a.index(row, col));
}
}
}
// restore the type
//
return this_a.changeType(old_type);
}
// explicit instantiations
//
#ifndef ISIP_TEMPLATE_complex
template boolean
MMatrixMethods::sub<ISIP_TEMPLATE_TARGET, Byte, byte>
(MMatrix<ISIP_TEMPLATE_TARGET>&, const MMatrix<ISIP_TEMPLATE_TARGET>&,
const MMatrix<Byte, byte>&);
template boolean
MMatrixMethods::sub<ISIP_TEMPLATE_TARGET, Ushort, ushort>
(MMatrix<ISIP_TEMPLATE_TARGET>&, const MMatrix<ISIP_TEMPLATE_TARGET>&,
const MMatrix<Ushort, ushort>&);
template boolean
MMatrixMethods::sub<ISIP_TEMPLATE_TARGET, Ulong, ulong>
(MMatrix<ISIP_TEMPLATE_TARGET>&, const MMatrix<ISIP_TEMPLATE_TARGET>&,
const MMatrix<Ulong, ulong>&);
template boolean
MMatrixMethods::sub<ISIP_TEMPLATE_TARGET, Ullong, ullong>
(MMatrix<ISIP_TEMPLATE_TARGET>&, const MMatrix<ISIP_TEMPLATE_TARGET>&,
const MMatrix<Ullong, ullong>&);
template boolean
MMatrixMethods::sub<ISIP_TEMPLATE_TARGET, Short, short>
(MMatrix<ISIP_TEMPLATE_TARGET>&, const MMatrix<ISIP_TEMPLATE_TARGET>&,
const MMatrix<Short, short>&);
template boolean
MMatrixMethods::sub<ISIP_TEMPLATE_TARGET, Long, long>
(MMatrix<ISIP_TEMPLATE_TARGET>&, const MMatrix<ISIP_TEMPLATE_TARGET>&,
const MMatrix<Long, long>&);
template boolean
MMatrixMethods::sub<ISIP_TEMPLATE_TARGET, Llong, llong>
(MMatrix<ISIP_TEMPLATE_TARGET>&, const MMatrix<ISIP_TEMPLATE_TARGET>&,
const MMatrix<Llong, llong>&);
template boolean
MMatrixMethods::sub<ISIP_TEMPLATE_TARGET, Float, float>
(MMatrix<ISIP_TEMPLATE_TARGET>&, const MMatrix<ISIP_TEMPLATE_TARGET>&,
const MMatrix<Float, float>&);
template boolean
MMatrixMethods::sub<ISIP_TEMPLATE_TARGET, Double, double>
(MMatrix<ISIP_TEMPLATE_TARGET>&, const MMatrix<ISIP_TEMPLATE_TARGET>&,
const MMatrix<Double, double>&);
#endif
#ifdef ISIP_TEMPLATE_complex
template boolean
MMatrixMethods::sub<ISIP_TEMPLATE_TARGET, ComplexLong, complexlong>
(MMatrix<ISIP_TEMPLATE_TARGET>&, const MMatrix<ISIP_TEMPLATE_TARGET>&,
const MMatrix<ComplexLong, complexlong>&);
template boolean
MMatrixMethods::sub<ISIP_TEMPLATE_TARGET, ComplexFloat, complexfloat>
(MMatrix<ISIP_TEMPLATE_TARGET>&, const MMatrix<ISIP_TEMPLATE_TARGET>&,
const MMatrix<ComplexFloat, complexfloat>&);
template boolean
MMatrixMethods::sub<ISIP_TEMPLATE_TARGET, ComplexDouble, complexdouble>
(MMatrix<ISIP_TEMPLATE_TARGET>&, const MMatrix<ISIP_TEMPLATE_TARGET>&,
const MMatrix<ComplexDouble, complexdouble>&);
#endif
[snipped]
//
// end of file
Up |
Home |
Courses |
Projects |
Proposals |
Publications
Please direct questions or comments to
joseph.picone@isip.piconepress.com