cmatrix.c

大地测量专业计算软件

          }
          return TRUE;
        }
        else
        {
          return TRUE;
        }
      }
      else if( M->nrows == nrows && M->ncols < ncols )
      {
        if( setToZero )
        {
          if( !MTX_Zero( M ) )
          {
            MTX_ERROR_MSG( "MTX_Zero returned FALSE." );
            return FALSE;
          }
        }
        if( !MTX_AddZeroValuedColumns( M, ncols-M->ncols ) )
        {
          MTX_ERROR_MSG( "MTX_Zero returned FALSE." );
          return FALSE;
        }
        return TRUE;
      }
      else if( M->nrows == nrows && M->ncols > ncols )
      {
        if( MTX_RemoveColumnsAfterIndex( M, ncols-1 ) == FALSE )
        {
          MTX_ERROR_MSG( "MTX_RemoveColumnsAfterIndex returned FALSE." );
          return FALSE;
        }
        if( setToZero )
        {
          if( !MTX_Zero( M ) )
          {
            MTX_ERROR_MSG( "MTX_Zero returned FALSE." );
            return FALSE;
          }
        }
        return TRUE;
      }
    }    
  }

  // The matrix must be built from scratch.
  MTX_Free( M );

  M->isReal = isReal;
  M->nrows = nrows;
  M->ncols = 0;

  // allocate the column array
  if( isReal )
  {
    M->data = (double**)malloc( ncols*sizeof(double*) );
    if( !M->data )
    {
      MTX_ERROR_MSG( "malloc returned NULL." );
      return FALSE;
    }
  }
  else
  {
    M->cplx = (stComplex**)malloc( ncols*sizeof(stComplex*) );
    if( !M->cplx )
    {
      MTX_ERROR_MSG( "malloc returned NULL." );
      return FALSE;
    }
  }


  // for each column allocate the rows
  if( isReal )
  {
    for( j = 0; j < ncols; j++ )
    {
      if( setToZero )
        M->data[j] = (double*)calloc( nrows, sizeof(double) );
      else
        M->data[j] = (double*)malloc( nrows*sizeof(double) );
      if( !M->data[j] )
      {
        // this is most likely to occur if allocating more memory than available
        MTX_ERROR_MSG( "malloc or calloc returned NULL." );
        MTX_Free( M );
        return FALSE;
      }
      M->ncols++;
    }
  }
  else
  {
    for( j = 0; j < ncols; j++ )
    {
      M->cplx[j] = (stComplex*)malloc( nrows*sizeof(stComplex) );
      if( !M->cplx[j] )
      {
        // this is most likely to occur if allocating more memory than available
        MTX_ERROR_MSG( "malloc returned NULL." );
        MTX_Free( M );
        return FALSE;
      }
      if( setToZero )
      {
        for( i = 0; i < nrows; i++ )
        {
          M->cplx[j][i].re = 0.0;
          M->cplx[j][i].im = 0.0;
        }
      }
      M->ncols++;
    }    
  }

  return TRUE;
}

// Set a scalar value in the matrix.
BOOL MTX_SetValue( MTX *M, const unsigned row, const unsigned col, const double value )
{
  if( MTX_isNull( M ) )
  {
    MTX_ERROR_MSG( "NULL Matrix" );
    return FALSE;
  }

  if( row >= M->nrows )
  {
    MTX_ERROR_MSG( "if( row >= M->nrows )" );
    return FALSE;
  }

  if( col >= M->ncols )
  {
    MTX_ERROR_MSG( "if( col >= M->ncols )" );
    return FALSE;
  }

  if( M->isReal )
  {
    M->data[col][row] = value;
  }
  else
  {
    M->cplx[col][row].re = value;
    M->cplx[col][row].im = 0.0;
  }

  return TRUE;
}

// Set a complex value in the matrix.
BOOL MTX_SetComplexValue( MTX *M, const unsigned row, const unsigned col, const double re, const double im )
{
  if( MTX_isNull( M ) )
  {
    MTX_ERROR_MSG( "NULL Matrix" );
    return FALSE;
  }

  if( row >= M->nrows )
  {
    MTX_ERROR_MSG( "if( row >= M->nrows )" );
    return FALSE;
  }

  if( col >= M->ncols )
  {
    MTX_ERROR_MSG( "if( col >= M->ncols )" );
    return FALSE;
  }

  if( M->isReal && im != 0.0 )
  {
    if( !MTX_ConvertRealToComplex( M ) )
    {
      MTX_ERROR_MSG( "MTX_ConvertRealToComplex returned FALSE." );
      return FALSE;
    }
  }

  if( M->isReal )
  {
    // im == 0.0
    M->data[col][row] = re;
  }
  else
  {
    M->cplx[col][row].re = re;
    M->cplx[col][row].im = im;
  }

  return TRUE;
}


// Matrix M = Re + Im*i, where Re and Im are real matrices.
BOOL MTX_Complex( MTX *M, const MTX *Re, const MTX *Im )
{
  unsigned i = 0;
  unsigned j = 0;

  if( MTX_isNull( Re ) )
  {
    MTX_ERROR_MSG( "NULL Matrix" );
    return FALSE;
  }
  if( MTX_isNull( Im ) )
  {
    MTX_ERROR_MSG( "NULL Matrix" );
    return FALSE;
  }
  if( !MTX_isSameSize( Re, Im ) )
  {
    MTX_ERROR_MSG( "MTX_isSameSize returned FALSE." );
    return FALSE;
  }

  if( !Re->isReal )
  {
    MTX_ERROR_MSG( "if( !Re->isReal )" );
    return FALSE;
  }
  if( !Im->isReal )
  {
    MTX_ERROR_MSG( "if( !Im->isReal )" );
    return FALSE;
  }

  if( M->isReal )
  {
    if( !MTX_Malloc( M, Re->nrows, Re->ncols, FALSE ) )
    {
      MTX_ERROR_MSG( "MTX_Malloc returned FALSE." );
      return FALSE;
    }
  }
  else
  {
    if( !MTX_Resize( M, Re->nrows, Re->ncols, FALSE ) )
    {
      MTX_ERROR_MSG( "MTX_Resize returned FALSE." );
      return FALSE;
    }
  }

  for( j = 0; j < M->ncols; j++ )
  {
    for( i = 0; i < M->nrows; i++ )
    {
      M->cplx[j][i].re = Re->data[j][i];
      M->cplx[j][i].im = Im->data[j][i];
    }
  }
  return TRUE;
}


// Set the specified column in Matrix M to Re + Im*i, where Re and Im are real matrices.
// The dimensions of M must already be valid.
BOOL MTX_SetComplexColumn( MTX *M, const unsigned col, const MTX *Re, const MTX *Im )
{
  unsigned i = 0;
  unsigned j = 0;

  // check that M is complex
  if( M->isReal )
  {
    MTX_ERROR_MSG( "if( M->isReal )" );
    return FALSE;
  }

  if( MTX_isNull( Re ) )
  {
    MTX_ERROR_MSG( "NULL Matrix" );
    return FALSE;
  }
  if( MTX_isNull( Im ) )
  {
    MTX_ERROR_MSG( "NULL Matrix" );
    return FALSE;
  }
  if( Re->ncols != 1 )
  {
    MTX_ERROR_MSG( "if( Re->ncols != 1 )" );
    return FALSE;
  }
  if( !MTX_isSameSize( Re, Im ) )
  {
    MTX_ERROR_MSG( "MTX_isSameSize returned FALSE." );
    return FALSE;
  }

  if( !Re->isReal )
  {
    MTX_ERROR_MSG( "if( !Re->isReal )" );
    return FALSE;
  }
  if( !Im->isReal )
  {
    MTX_ERROR_MSG( "if( !Im->isReal )" );
    return FALSE;
  }

  // check that M has the right dimension
  if( M->nrows != Re->nrows )
  {
    MTX_ERROR_MSG( "if( M->nrows != Re->nrows )" );
    return FALSE;
  }
  if( col >= M->ncols )
  {
    MTX_ERROR_MSG( "if( col >= M->ncols )" );
    return FALSE;
  }

  for( i = 0; i < M->nrows; i++ )
  {
    M->cplx[col][i].re = Re->data[0][i];
    M->cplx[col][i].im = Im->data[0][i];
  }
  return TRUE;
}


BOOL MTX_ConvertRealToComplex( MTX *M )
{
  unsigned i = 0;
  unsigned j = 0;
  unsigned k = 0;

  if( MTX_isNull( M ) )
  {
    MTX_ERROR_MSG( "NULL Matrix" );
    return FALSE;
  }

  if( !M->isReal )
    return TRUE; // already complex, nothing to do

  // allocate the complex column vector pointers
  M->cplx = (stComplex**)malloc( M->ncols*sizeof(stComplex*) );
  if( !M->cplx )
  {
    MTX_ERROR_MSG( "malloc retuned NULL." );
    return FALSE;
  }

  for( j = 0; j < M->ncols; j++ )
  {
    M->cplx[j] = (stComplex*)malloc( sizeof(stComplex)*(M->nrows) );
    if( !(M->cplx[j]) )
    {
      // this is most likely to occur if allocating more memory than available
      for( k = 0; k < j; k++ ) // delete the complex column data already allocated
      {
        free( M->cplx[k] );
      }
      // delete the complex column pointer array
      free( M->cplx );
      M->cplx = NULL;
      MTX_ERROR_MSG( "malloc retuned NULL." );
      return FALSE; // note, the matrix M is still valid as a real matrix
    }
    // now copy the real data to the complex
    for( i = 0; i < M->nrows; i++ )
    {
      M->cplx[j][i].re = M->data[j][i];
      M->cplx[j][i].im = 0.0;
    }
  }

  // free the real data
  for( j = 0; j < M->ncols; j++ )
  {
    free( M->data[j] );
  }
  // free the array of real pointers
  free( M->data );
  M->data = NULL;
  M->isReal = FALSE;

  // successfully converted the matrix from real to complex
  return TRUE;
}


BOOL MTX_ConvertComplexToReal( MTX *M )
{
  return MTX_static_ConvertComplexTo( M, TRUE );
}

BOOL MTX_ConvertComplexToImag( MTX *M )
{
  return MTX_static_ConvertComplexTo( M, FALSE );
}

BOOL MTX_static_ConvertComplexTo( MTX *M, BOOL useReal )
{
  unsigned i = 0;
  unsigned j = 0;
  unsigned k = 0;

  if( MTX_isNull( M ) )
  {
    MTX_ERROR_MSG( "NULL Matrix" );
    return FALSE;
  }

  // Deal with special case of already real values.
  if( useReal )
  {
    if( M->isReal )
      return TRUE; // already real, nothing to do
  }

  // Deal with special case of trying to get complex values from a real matrix.
  if( !useReal && M->isReal )
  {
    if( !MTX_Zero(M) )
    {
      MTX_ERROR_MSG( "MTX_Zero returned FALSE." );
      return FALSE;
    }
    return TRUE;
  }

  // allocate the complex column vector pointers
  M->data = (double**)malloc( (M->ncols)*sizeof(double*) );
  if( !M->data )
  {
    MTX_ERROR_MSG( "if( !M->data )" );
    return FALSE;
  }

  for( j = 0; j < M->ncols; j++ )
  {
    M->data[j] = (double*)malloc( sizeof(double)*(M->nrows) );
    if( !(M->data[j]) )
    {
      // this is most likely to occur if allocating more memory than available
      for( k = 0; k < j; k++ ) // delete the real column data already allocated
      {
        free( M->data[k] );
      }
      // delete the real column pointer array
      free( M->data );
      M->data = NULL;
      MTX_ERROR_MSG( "malloc returned NULL." );
      return FALSE; // note, the matrix M is still valid as a complex matrix
    }
    // now copy the complex real component of the data to the real data.
    for( i = 0; i < M->nrows; i++ )
    {
      if( useReal )
        M->data[j][i] = M->cplx[j][i].re;
      else
        M->data[j][i] = M->cplx[j][i].im;
    }
  }

  // free the complex data
  for( j = 0; j < M->ncols; j++ )
  {
    free( M->cplx[j] );
  }
  // free the array of real pointers
  free( M->cplx );
  M->cplx = NULL;
  M->isReal = TRUE;

  // successfully converted the matrix from complex to real
  return TRUE;
}


// Extract the real component of matrix M
BOOL MTX_Real( const MTX *M, MTX *Re )
{
  unsigned i = 0;
  unsigned j = 0;

  if( MTX_isNull( M ) )
  {
    MTX_ERROR_MSG( "NULL Matrix" );
    return FALSE;
  }

  if( M->isReal )
    return MTX_Copy( M, Re );

  if( !MTX_Malloc( Re, M->nrows, M->ncols, TRUE ) )
  {
    MTX_ERROR_MSG( "MTX_Malloc returned FALSE." );
    return FALSE;
  }

  for( j = 0; j < M->ncols; j++ )
  {
    for( i = 0; i < M->nrows; i++ )
    {
      Re->data[j][i] = M->cplx[j][i].re;
    }
  }
  return TRUE;
}


BOOL MTX_isReal( MTX *M, BOOL *isReal )
{
  if( MTX_isNull( M ) )
  {
    *isReal = TRUE;
    return TRUE; // A null matrix is real by default.
  }

  if( M->isReal )
  {
    *isReal = TRUE;
  }
  else
  {
    double maxabs;
    MTX imagM;

    MTX_Init(&imagM);

    if( !MTX_Imag(M,&imagM) )
    {
      MTX_ERROR_MSG( "MTX_Imag returned FALSE." );
      return FALSE;
    }

    if( !MTX_MaxAbs(&imagM,&maxabs) )
    {
      MTX_ERROR_MSG( "MTX_MaxAbs returned FALSE." );
      return FALSE;
    }

    if( maxabs == 0.0 )
    {
      if( !MTX_ConvertComplexToReal(M) )
      {
        MTX_ERROR_MSG( "MTX_ConvertComplexToReal returned FALSE." );
        return FALSE;
      }
      *isReal = TRUE;
    }
    else
    {
      *isReal = FALSE;
    }

    MTX_Free(&imagM);
  }

  return TRUE;
}

BOOL MTX_RealColumn( const MTX *M, const unsigned col, MTX *Re )
{
  unsigned i = 0;

  if( MTX_isNull( M ) )
  {
    MTX_ERROR_MSG( "NULL Matrix" );
    return FALSE;
  }

  if( col >= M->ncols )
  {
    MTX_ERROR_MSG( "if( col >= M->ncols )" );
    return FALSE;
  }

  if( M->isReal )
    return MTX_CopyColumn( M, col, Re );

  if( !MTX_Malloc( Re, M->nrows, 1, TRUE ) )
  {