// file: $isip/class/algo/Math/math_08.cc
// version: $Id: math_08.cc 8037 2002-04-01 20:04:33Z gao $
//

// isip include files
//
#include "Math.h"

// method: computeFuncCalcEnumerateDouble
//
// arguments:
//  AlgorithmData& output: (output) output data
//  const Vector<AlgorithmData>& input: (input) input data
//
// return: a bool8 value indicating status
//
// this method implements following function:
//   Y = operation(0) {a(0)*function(0)(X(0))} operation(1)
//   {a(1)*function(1)(X(1))} operation(2) {a(2)*function(2)(X(2))}...
//   ... + const,
// where X(0), X(1), ... are the parts of input VectorDouble splitted
// by number of operands. These can be VectorDouble. If
// user wants to implement this for scalar or single vector, then
// number of operands will be 1. For scalar input the length of the
// input vector will be 1.
//
bool8 Math::computeFuncCalcEnumerateDouble(AlgorithmData& output_a,
					     const Vector<AlgorithmData>&
					     input_a) {

  // check algorithm and implementation names
  //
  if ((algorithm_d != FUNCTION_CALCULATOR) ||
      (implementation_d != ENUMERATE)) {
    return Error::handle(name(), L"computeFuncCalcEnumerateDouble", Error::ARG,
                         __FILE__, __LINE__);
  }

  // compute the number of combinations
  //
  if (num_operands_d != input_a.length()) {
    return Error::handle(name(), L"computeFuncCalcEnumerateDouble",
			 Error::ARG, __FILE__, __LINE__);
  }

  // input should be either vector or matrix
  //
  for (int32 j = 0; j < num_operands_d; j++) {
    if ((input_a(j).getDataType() != AlgorithmData::VECTOR_DOUBLE) &&
	(input_a(j).getDataType() != AlgorithmData::MATRIX_DOUBLE)) {
      return Error::handle(name(), L"computeFuncCalcEnumerateDouble",
			   Error::ARG, __FILE__, __LINE__);
    }
  }

  // set the mode of output
  //
  if (input_a(0).getDataType() == AlgorithmData::VECTOR_DOUBLE) {
    output_a.makeVectorDouble();
  }
  else {
    output_a.makeMatrixDouble();
  }

  // loop over the number of input operands
  //
  for (int32 j = 0; j < num_operands_d; j++) {
    
    // declare temporary input variable
    //
    AlgorithmData tmp_input;
    
    // set the type to either matrix or vector, depending upon the
    // type of the input
    //
    if (input_a(j).getDataType() == AlgorithmData::VECTOR_DOUBLE) {
      tmp_input.makeVectorDouble().assign(input_a(j).getVectorDouble());
    }
    else {
      tmp_input.makeMatrixDouble().assign(input_a(j).getMatrixDouble());
    }

    // branch on function:
    //  apply function and corresponding weight
    //
    if (function_d(j) == (int32)IDENTITY) {
      
      // do nothing
      //
    }
    
    // function: EXP
    //
    else if (function_d(j) == (int32)EXP) {
      if (input_a(j).getDataType() != AlgorithmData::VECTOR_DOUBLE) {
	return Error::handle(name(), L"computeFuncCalcEnumerateDouble",
			     ERR_FNDTYP, __FILE__, __LINE__);
      }
      tmp_input.getVectorDouble().exp();
    }
    
    // function: EXP2
    //    
    else if (function_d(j) == (int32)EXP2) {
      if (input_a(j).getDataType() != AlgorithmData::VECTOR_DOUBLE) {
	return Error::handle(name(), L"computeFuncCalcEnumerateDouble",
			     ERR_FNDTYP, __FILE__, __LINE__);
      }
      tmp_input.getVectorDouble().exp2();
    }
    
    // function: EXP10
    //    
    else if (function_d(j) == (int32)EXP10) {
      if (input_a(j).getDataType() != AlgorithmData::VECTOR_DOUBLE) {
	return Error::handle(name(), L"computeFuncCalcEnumerateDouble",
			     ERR_FNDTYP, __FILE__, __LINE__);
      }
      tmp_input.getVectorDouble().exp10();
    }
    
    // function: FACTORIAL
    //    
    else if (function_d(j) == (int32)FACTORIAL) {
      if (input_a(j).getDataType() != AlgorithmData::VECTOR_DOUBLE) {
	return Error::handle(name(), L"computeFuncCalcEnumerateDouble",
			     ERR_FNDTYP, __FILE__, __LINE__);
      }
      tmp_input.getVectorDouble().factorial();
    }
    
    // function: LOG
    //        
    else if (function_d(j) == (int32)LOG) {
      if (input_a(j).getDataType() != AlgorithmData::VECTOR_DOUBLE) {
	return Error::handle(name(), L"computeFuncCalcEnumerateDouble",
			     ERR_FNDTYP, __FILE__, __LINE__);
      }
  
      // check vector elements
      //
      for (int32 index =  tmp_input.getVectorDouble().length() - 1;
	   index >= 0; index--) {
	if (tmp_input.getVectorDouble()(index) <= (float64)0.0 ) {
	  return Error::handle(name(), L"computeFuncCalcEnumerateDouble",
			       Error::ARG, __FILE__, __LINE__);
	}
      }
      tmp_input.getVectorDouble().log();
    }
    
    // function: LOG2
    //            
    else if (function_d(j) == (int32)LOG2) {
      if (input_a(j).getDataType() != AlgorithmData::VECTOR_DOUBLE) {
	return Error::handle(name(), L"computeFuncCalcEnumerateDouble",
			     ERR_FNDTYP, __FILE__, __LINE__);
      }
      
      // check vector elements
      //
      for (int32 index =  tmp_input.getVectorDouble().length() - 1;
	   index >= 0; index--) {
	if (tmp_input.getVectorDouble()(index) <= (float64)0.0 ) {
	  return Error::handle(name(), L"computeFuncCalcEnumerateDouble",
			       Error::ARG, __FILE__, __LINE__);
	}
      }
      tmp_input.getVectorDouble().log2();
    }
    
    // function: LOG10
    //            
    else if (function_d(j) == (int32)LOG10) {
      if (input_a(j).getDataType() != AlgorithmData::VECTOR_DOUBLE) {
	return Error::handle(name(), L"computeFuncCalcEnumerateDouble",
			     ERR_FNDTYP, __FILE__, __LINE__);
      }
      
      // check vector elements
      //
      for (int32 index =  tmp_input.getVectorDouble().length() - 1;
	   index >= 0; index--) {
	if (tmp_input.getVectorDouble()(index) <= (float64)0.0 ) {
	  return Error::handle(name(), L"computeFuncCalcEnumerateDouble",
			       Error::ARG, __FILE__, __LINE__);
	}
      }
      tmp_input.getVectorDouble().log10();
    }
    
    // function: LOG1P
    //            
    else if (function_d(j) == (int32)LOG1P) {
      if (input_a(j).getDataType() != AlgorithmData::VECTOR_DOUBLE) {
	return Error::handle(name(), L"computeFuncCalcEnumerateDouble",
			     ERR_FNDTYP, __FILE__, __LINE__);
      }
      
      // check vector elements
      //
      for (int32 index =  tmp_input.getVectorDouble().length() - 1;
	   index >= 0; index--) {
	if (tmp_input.getVectorDouble()(index) <= (float64)0.0 ) {
	  return Error::handle(name(), L"computeFuncCalcEnumerateDouble",
			       Error::ARG, __FILE__, __LINE__);
	}
      }
      tmp_input.getVectorDouble().log1p();
    }
    
    // function: ABS
    //            
    else if (function_d(j) == (int32)ABS) {
      if (input_a(j).getDataType() != AlgorithmData::VECTOR_DOUBLE) {
	return Error::handle(name(), L"computeFuncCalcEnumerateDouble",
			     ERR_FNDTYP, __FILE__, __LINE__);
      }
      tmp_input.getVectorDouble().abs();
    }
    
    // function: NEG
    //            
    else if (function_d(j) == (int32)NEG) {
      if (input_a(j).getDataType() != AlgorithmData::VECTOR_DOUBLE) {
	tmp_input.getVectorDouble().neg();
      }
      else if (input_a(j).getDataType() != AlgorithmData::MATRIX_DOUBLE) {
	tmp_input.getMatrixDouble().neg();
      }
      else {
	return Error::handle(name(), L"computeFuncCalcEnumerateDouble",
			     Error::ARG, __FILE__, __LINE__);
      }
    }
    
    // function: ROUND
    //                 
    else if (function_d(j) == (int32)ROUND) {
      if (input_a(j).getDataType() != AlgorithmData::VECTOR_DOUBLE) {
	return Error::handle(name(), L"computeFuncCalcEnumerateDouble",
			     ERR_FNDTYP, __FILE__, __LINE__);
      }
      tmp_input.getVectorDouble().round();
    }
    
    // function: CEIL
    //                
    else if (function_d(j) == (int32)CEIL) {
      if (input_a(j).getDataType() != AlgorithmData::VECTOR_DOUBLE) {
	return Error::handle(name(), L"computeFuncCalcEnumerateDouble",
			     ERR_FNDTYP, __FILE__, __LINE__);
      }
      tmp_input.getVectorDouble().ceil();
    }
    
    // function: FLOOR
    //                    
    else if (function_d(j) == (int32)FLOOR) {
      if (input_a(j).getDataType() != AlgorithmData::VECTOR_DOUBLE) {
	return Error::handle(name(), L"computeFuncCalcEnumerateDouble",
			     ERR_FNDTYP, __FILE__, __LINE__);
      }
      tmp_input.getVectorDouble().floor();
    }
    
    // function: RFLOOR
    //                        
    else if (function_d(j) == (int32)RFLOOR) {
      if (input_a(j).getDataType() != AlgorithmData::VECTOR_DOUBLE) {
	return Error::handle(name(), L"computeFuncCalcEnumerateDouble",
			     ERR_FNDTYP, __FILE__, __LINE__);
      }
      tmp_input.getVectorDouble().rfloor();
    }
    
    // function: SIN
    //                        
    else if (function_d(j) == (int32)SIN) {
      if (input_a(j).getDataType() != AlgorithmData::VECTOR_DOUBLE) {
	return Error::handle(name(), L"computeFuncCalcEnumerateDouble",
			     ERR_FNDTYP, __FILE__, __LINE__);
      }
      tmp_input.getVectorDouble().sin();
    }
    
    // function: COS
    //                            
    else if (function_d(j) == (int32)COS) {
      if (input_a(j).getDataType() != AlgorithmData::VECTOR_DOUBLE) {
	return Error::handle(name(), L"computeFuncCalcEnumerateDouble",
			     ERR_FNDTYP, __FILE__, __LINE__);
      }
      tmp_input.getVectorDouble().cos();
    }

    // function: TAN
    //                            
    else if (function_d(j) == (int32)TAN) {
      if (input_a(j).getDataType() != AlgorithmData::VECTOR_DOUBLE) {
	return Error::handle(name(), L"computeFuncCalcEnumerateDouble",
			     ERR_FNDTYP, __FILE__, __LINE__);
      }
      tmp_input.getVectorDouble().tan();
    }
    
    // function: ASIN
    //                            
    else if (function_d(j) == (int32)ASIN) {
      if (input_a(j).getDataType() != AlgorithmData::VECTOR_DOUBLE) {
	return Error::handle(name(), L"computeFuncCalcEnumerateDouble",
			     ERR_FNDTYP, __FILE__, __LINE__);
      }
      tmp_input.getVectorDouble().asin();
    }
    
    // function: ACOS
    //                            
    else if (function_d(j) == (int32)ACOS) {
      if (input_a(j).getDataType() != AlgorithmData::VECTOR_DOUBLE) {
	return Error::handle(name(), L"computeFuncCalcEnumerateDouble",
			     ERR_FNDTYP, __FILE__, __LINE__);
      }
      tmp_input.getVectorDouble().acos();
    }
    
    // function: ATAN
    //                            
    else if (function_d(j) == (int32)ATAN) {
      if (input_a(j).getDataType() != AlgorithmData::VECTOR_DOUBLE) {
	return Error::handle(name(), L"computeFuncCalcEnumerateDouble",
			     ERR_FNDTYP, __FILE__, __LINE__);
      }
      tmp_input.getVectorDouble().atan();
    }
    
    // function: SQUARE
    //                            
    else if (function_d(j) == (int32)SQUARE) {
      if (input_a(j).getDataType() != AlgorithmData::VECTOR_DOUBLE) {
	return Error::handle(name(), L"computeFuncCalcEnumerateDouble",
			     ERR_FNDTYP, __FILE__, __LINE__);
      }
      tmp_input.getVectorDouble().square();
    }
    
    // function: SQRT
    //                            
    else if (function_d(j) == (int32)SQRT) {
      if (input_a(j).getDataType() != AlgorithmData::VECTOR_DOUBLE) {
	return Error::handle(name(), L"computeFuncCalcEnumerateDouble",
			     ERR_FNDTYP, __FILE__, __LINE__);
      }
      tmp_input.getVectorDouble().sqrt();
    }
    
    // function: INVERSE
    //                            
    else if (function_d(j) == (int32)INVERSE) {
      if (input_a(j).getDataType() != AlgorithmData::MATRIX_DOUBLE) {
	return Error::handle(name(), L"computeFuncCalcEnumerateDouble",
			     ERR_FNDTYP, __FILE__, __LINE__);
      }
      tmp_input.getMatrixDouble().inverse();
    }
    
    // function: TRANSPOSE
    //                            
    else if (function_d(j) == (int32)TRANSPOSE) {
      if (input_a(j).getDataType() != AlgorithmData::MATRIX_DOUBLE) {
	return Error::handle(name(), L"computeFuncCalcEnumerateDouble",
			     ERR_FNDTYP, __FILE__, __LINE__);
      }
      tmp_input.getMatrixDouble().transpose();
    }
    
    // invalid function
    //
    else {
      return Error::handle(name(), L"computeFuncCalcEnumerateDouble",
			   ERR_UNKFUN, __FILE__, __LINE__);
    }
    
    // apply weight
    //
    if (input_a(j).getDataType() == AlgorithmData::VECTOR_DOUBLE) {
      tmp_input.getVectorDouble().mult(weight_d(j));
    }
    else if (input_a(j).getDataType() == AlgorithmData::MATRIX_DOUBLE) {
      tmp_input.getMatrixDouble().mult(weight_d(j));
    }
    else {
      return Error::handle(name(), L"computeFuncCalcEnumerateDouble",
			   Error::ARG, __FILE__, __LINE__);
    }

    // branch on operation:
    //  Operation: ASSIGN
    //
    if (operation_d(j) == (int32)ASSIGN) {
      
      // this operation is possible between same data types
      //
      if (output_a.getDataType() != input_a(j).getDataType()) {
	return Error::handle(name(), L"computeFuncCalcEnumerateDouble",
			     ERR_MATCH, __FILE__, __LINE__);
      }
      if (output_a.getDataType() == AlgorithmData::VECTOR_DOUBLE) {
	output_a.makeVectorDouble().assign(tmp_input.getVectorDouble());
      }
      else {
	output_a.makeMatrixDouble().assign(tmp_input.getMatrixDouble());
      }
    }
    
    //  Operation: ADD
    //
    else if (operation_d(j) == (int32)ADD) {
      
      computeAddDouble(output_a, tmp_input);
    }

    //  Operation: SUBTRACT
    //
    else if (operation_d(j) == (int32)SUBTRACT) {
      
      computeSubDouble(output_a, tmp_input);
    }

    //  Operation: MULTIPLY
    //
    else if (operation_d(j) == (int32)MULTIPLY) {

      computeMultDouble(output_a, tmp_input);
    }
    
    //  Operation: DIVIDE
    //
    else if (operation_d(j) == (int32)DIVIDE) {
      
      computeDivDouble(output_a, tmp_input);
    }
    
    // invalid operation
    //
    else {
      return Error::handle(name(), L"computeFuncCalcEnumerateDouble",
			   ERR_UNKOPE, __FILE__, __LINE__);
    }
  }

  // add the constant term
  //
  if (output_a.getDataType() == AlgorithmData::VECTOR_DOUBLE) {
    output_a.getVectorDouble().add((float64)const_d);
  }
  else {
    output_a.getMatrixDouble().add((float64)const_d);
  }
  
  // exit gracefully
  //
  return true;
}

// method: computeAddDouble
//
// arguments:
//  AlgorithmData& output: (input/output) 1st operand and output
//  const AlgorithmData& input: (input) 2nd operand
//
// return: logical error status
//
bool8 Math::computeAddDouble(AlgorithmData& output_a,
			       const AlgorithmData& input_a) const {
  
  if (input_a.getDataType() == AlgorithmData::VECTOR_DOUBLE) {
    
    int32 ilen = input_a.getVectorDouble().length();

    if (output_a.getDataType() == AlgorithmData::VECTOR_DOUBLE) {

      int32 olen = output_a.getVectorDouble().length();
      
      // if either input is of length 1, perform scalar addition
      //
      if (ilen == 1) {
	float64 scalar = input_a.getVectorDouble()(0);
	return output_a.getVectorDouble().add(scalar);
      }
      else if (olen == 1) {
	float64 scalar = output_a.getVectorDouble()(0);
	return output_a.getVectorDouble().add(input_a.getVectorDouble(), scalar);
      }
      
      // if they are the same length, do vector addition
      //
      else if (ilen == olen) {
	return output_a.getVectorDouble().add(input_a.getVectorDouble());
      }

      // else error
      //
      else {
	return Error::handle(name(), L"computeAddDouble", ERR_MATCH,
			     __FILE__, __LINE__);
      }
    }

    // 1 vector, 1 matrix. if the vector is of length 1, treat it as a
    // scalar
    //
    else if (output_a.getDataType() == AlgorithmData::MATRIX_DOUBLE) {

      if (ilen == 1) {
	float64 scalar = input_a.getVectorDouble()(0);
	return output_a.getMatrixDouble().add(scalar);
      }
      else {
	return Error::handle(name(), L"computeAddDouble", ERR_MATCH,
			     __FILE__, __LINE__);
      }
    }
  }

  else if (input_a.getDataType() == AlgorithmData::MATRIX_DOUBLE) {

    // 1 vector, 1 matrix. if the vector is of length 1, treat it as a
    // scalar
    //
    if (output_a.getDataType() == AlgorithmData::VECTOR_DOUBLE) {
      int32 olen = output_a.getVectorDouble().length();
      if (olen == 1) {
	float64 scalar = output_a.getVectorDouble()(0);
	return output_a.makeMatrixDouble().add(input_a.getMatrixDouble(), scalar);
      }
      else {
	return Error::handle(name(), L"computeAddDouble", ERR_MATCH, __FILE__,
			     __LINE__);
      }
    }

    // if both are matrices, just add them
    //
    else if (output_a.getDataType() == AlgorithmData::MATRIX_DOUBLE) {
      return output_a.getMatrixDouble().add(input_a.getMatrixDouble());
    }
  }

  // exit gracefully
  //
  return true;
}


// method: computeSubDouble
//
// arguments:
//  AlgorithmData& output: (input/output) 1st operand and output
//  const AlgorithmData& input: (input) 2nd operand
//
// return: logical error status
//
bool8 Math::computeSubDouble(AlgorithmData& output_a,
			       const AlgorithmData& input_a) const {
  
  if (input_a.getDataType() == AlgorithmData::VECTOR_DOUBLE) {
    
    int32 ilen = input_a.getVectorDouble().length();

    if (output_a.getDataType() == AlgorithmData::VECTOR_DOUBLE) {

      int32 olen = output_a.getVectorDouble().length();
      
      // if either input is of length 1, perform scalar subtraction
      //
      if (ilen == 1) {
	float64 scalar = input_a.getVectorDouble()(0);
	return output_a.getVectorDouble().sub(scalar);
      }
      else if (olen == 1) {

	// note this is different from addition since the operation is
	// not cummative
	//
	float64 scalar = output_a.getVectorDouble()(0);
	output_a.getVectorDouble().setLength(ilen);
	output_a.getVectorDouble().assign(scalar);
	return output_a.getVectorDouble().sub(input_a.getVectorDouble());
      }
      
      // if they are the same length, do vector subtraction
      //
      else if (ilen == olen) {
	return output_a.getVectorDouble().sub(input_a.getVectorDouble());
      }

      // else error
      //
      else {
	return Error::handle(name(), L"computeSubDouble", ERR_MATCH,
			     __FILE__, __LINE__);
      }
    }

    // 1 vector, 1 matrix. if the vector is of length 1, treat it as a
    // scalar
    //
    else if (output_a.getDataType() == AlgorithmData::MATRIX_DOUBLE) {

      if (ilen == 1) {
	float64 scalar = input_a.getVectorDouble()(0);
	return output_a.getMatrixDouble().sub(scalar);
      }
      else {
	return Error::handle(name(), L"computeSubDouble", ERR_MATCH,
			     __FILE__, __LINE__);
      }
    }
  }

  else if (input_a.getDataType() == AlgorithmData::MATRIX_DOUBLE) {

    // 1 vector, 1 matrix. if the vector is of length 1, treat it as a
    // scalar
    //
    if (output_a.getDataType() == AlgorithmData::VECTOR_DOUBLE) {
      int32 olen = output_a.getVectorDouble().length();
      if (olen == 1) {

	// note this is different from addition since the operation is
	// not cummative
	//
	float64 scalar = output_a.getVectorDouble()(0);
	output_a.makeMatrixDouble().setDimensions(input_a.getMatrixDouble());
	output_a.getMatrixDouble().assign(scalar);
	return output_a.getMatrixDouble().sub(input_a.getMatrixDouble());
      }
      else {
	return Error::handle(name(), L"computeSubDouble", ERR_MATCH, __FILE__,
			     __LINE__);
      }
    }

    // if both are matrices, just subtract them
    //
    else if (output_a.getDataType() == AlgorithmData::MATRIX_DOUBLE) {
      return output_a.getMatrixDouble().sub(input_a.getMatrixDouble());
    }
  }

  // exit gracefully
  //
  return true;
}

// method: computeMultDouble
//
// arguments:
//  AlgorithmData& output: (input/output) 1st operand and output
//  const AlgorithmData& input: (input) 2nd operand
//
// return: logical error status
//
bool8 Math::computeMultDouble(AlgorithmData& output_a,
				const AlgorithmData& input_a) const {
  
  if (input_a.getDataType() == AlgorithmData::VECTOR_DOUBLE) {
    
    int32 ilen = input_a.getVectorDouble().length();

    if (output_a.getDataType() == AlgorithmData::VECTOR_DOUBLE) {

      int32 olen = output_a.getVectorDouble().length();

      // if either input is of length 1, perform scalar multiplication
      //
      if (ilen == 1) {
	float64 scalar = input_a.getVectorDouble()(0);
	return output_a.getVectorDouble().mult(scalar);
      }
      else if (olen == 1) {

	float64 scalar = output_a.getVectorDouble()(0);
	return output_a.getVectorDouble().mult(input_a.getVectorDouble(), scalar);
      }
      
      // if they are the same length, do vector multiplication
      //
      else if (ilen == olen) {
	return output_a.getVectorDouble().mult(input_a.getVectorDouble());
      }

      // else error
      //
      else {
	return Error::handle(name(), L"computeMultDouble", ERR_MATCH,
			     __FILE__, __LINE__);
      }
    }
    else if (output_a.getDataType() == AlgorithmData::MATRIX_DOUBLE) {

      // if the vector is length 1, do scalar multiplication
      //
      if (ilen == 1) {
	float64 scalar = input_a.getVectorDouble()(0);
	return output_a.getMatrixDouble().mult(scalar);
      }

      // if the vector is the correct length for multv, do it
      //
      else if (ilen == output_a.getMatrixDouble().getNumColumns()) {
	VectorDouble tmp_out;
	output_a.getMatrixDouble().multv(tmp_out, input_a.getVectorDouble());
	output_a.makeVectorDouble().assign(tmp_out);
      }

      // else error
      //
      else {
	return Error::handle(name(), L"computeMultDouble", ERR_MATCH,
			     __FILE__, __LINE__);
      }
    }
  }

  // input is a matrix
  //
  else if (input_a.getDataType() == AlgorithmData::MATRIX_DOUBLE) {

    if (output_a.getDataType() == AlgorithmData::VECTOR_DOUBLE) {

      int32 olen = output_a.getVectorDouble().length();

      // if the vector is of length 1, do scalar multiplication
      //
      if (olen == 1) {
	float64 scalar = output_a.getVectorDouble()(0);
	output_a.makeMatrixDouble().assign(input_a.getMatrixDouble());
	return output_a.getMatrixDouble().mult(scalar);
      }

      // if the vector is the correct length for vmult, do it
      //
      else if (olen == input_a.getMatrixDouble().getNumRows()) {

	VectorDouble tmp_out;
	input_a.getMatrixDouble().vmult(tmp_out, output_a.getVectorDouble());
	output_a.getVectorDouble().assign(tmp_out);
      }
	
      // else error
      //
      else {
	return Error::handle(name(), L"computeMultDouble", ERR_MATCH,
			     __FILE__, __LINE__);
      }
    }
    else if (output_a.getDataType() == AlgorithmData::MATRIX_DOUBLE) {

      return output_a.getMatrixDouble().mult(input_a.getMatrixDouble());
    }
  }

  // exit gracefully
  //
  return true;
}

// method: computeDivDouble
//
// arguments:
//  AlgorithmData& output: (input/output) 1st operand and output
//  const AlgorithmData& input: (input) 2nd operand
//
// return: logical error status
//
bool8 Math::computeDivDouble(AlgorithmData& output_a,
			       const AlgorithmData& input_a) const {    

  
  // vector-matrix division is possible only if length of the
  // vector is 1. no other matrix division is possible
  //
  if ((input_a.getDataType() == AlgorithmData::VECTOR_DOUBLE) &&
      (output_a.getDataType() == AlgorithmData::MATRIX_DOUBLE)) {      

    if (input_a.getVectorDouble().length() > 1) {
      return Error::handle(name(), L"computeDivDouble",
			   ERR_MATCH, __FILE__, __LINE__);
    }
    VectorDouble vec(input_a.getVectorDouble());
    float64 scalar = vec(0);
    output_a.getMatrixDouble().div(scalar);
  }

  // both inputs are vectors
  //
  else if ((input_a.getDataType() == AlgorithmData::VECTOR_DOUBLE) &&
	   (output_a.getDataType() == AlgorithmData::VECTOR_DOUBLE)) {

    // if the input is of length 1, do scalar division
    //
    if (input_a.getVectorDouble().length() == 1) {
      float64 scalar = input_a.getVectorDouble()(0);

      return output_a.getVectorDouble().div(scalar);
    }

    // if the lengths are equal, do element-wise division
    //
    else if (input_a.getVectorDouble().length() == output_a.getVectorDouble().length()) {

      return output_a.getVectorDouble().div(input_a.getVectorDouble());

    } 
    else {
      return Error::handle(name(), L"computeDivDouble", ERR_MATCH,
			   __FILE__, __LINE__);
    }
  }
  else {
    return Error::handle(name(), L"computeDivDouble", ERR_MATCH, __FILE__,
			 __LINE__);
  }
  
  // exit gracefully
  //
  return true;
}