/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