matrix.c

spiht for linux this is used to decod and encode vedio i wich all enjoy

                 int num_rows, int num_cols)
{
  int return_value;
  QccVector col_vector1 = NULL;
  QccVector col_vector2 = NULL;
  int row, col;
  
  if ((input_block == NULL) ||
      (output_block == NULL))
    return(0);
  if ((num_rows <= 0) || (num_cols <= 0))
    return(0);
  
  if ((col_vector1 = QccVectorAlloc(num_rows)) == NULL)
    {
      QccErrorAddMessage("(QccMatrixDCT): Error calling QccVectorAlloc()");
      goto Error;
    }
  if ((col_vector2 = QccVectorAlloc(num_rows)) == NULL)
    {
      QccErrorAddMessage("(QccMatrixDCT): Error calling QccVectorAlloc()");
      goto Error;
    }
  
  for (row = 0; row < num_rows; row++)
    if (QccVectorDCT(input_block[row],
                     output_block[row],
                     num_cols))
      {
        QccErrorAddMessage("(QccMatrixDCT): Error calling QccVectorDCT()");
        goto Error;
      }
  
  for (col = 0; col < num_cols; col++)
    {
      for (row = 0; row < num_rows; row++)
        col_vector1[row] = output_block[row][col];
      if (QccVectorDCT(col_vector1,
                       col_vector2,
                       num_rows))
        {
          QccErrorAddMessage("(QccMatrixDCT): Error calling QccVectorDCT()");
          goto Error;
        }
      for (row = 0; row < num_rows; row++)
        output_block[row][col] = col_vector2[row];
    }
  
  return_value = 0;
  goto Return;
 Error:
  return_value = 1;
 Return:
  QccVectorFree(col_vector1);
  QccVectorFree(col_vector2);
  return(return_value);
}


int QccMatrixInverseDCT(const QccMatrix input_block,
                        QccMatrix output_block,
                        int num_rows, int num_cols)
{
  int return_value;
  QccVector col_vector1 = NULL;
  QccVector col_vector2 = NULL;
  int row, col;
  
  if ((input_block == NULL) ||
      (output_block == NULL))
    return(0);
  if ((num_rows <= 0) || (num_cols <= 0))
    return(0);
  
  if ((col_vector1 = QccVectorAlloc(num_rows)) == NULL)
    {
      QccErrorAddMessage("(QccMatrixInverseDCT): Error calling QccVectorAlloc()");
      goto Error;
    }
  if ((col_vector2 = QccVectorAlloc(num_rows)) == NULL)
    {
      QccErrorAddMessage("(QccMatrixInverseDCT): Error calling QccVectorAlloc()");
      goto Error;
    }
  
  for (row = 0; row < num_rows; row++)
    if (QccVectorInverseDCT(input_block[row],
                            output_block[row],
                            num_cols))
      {
        QccErrorAddMessage("(QccMatrixInverseDCT): Error calling QccVectorInverseDCT()");
        goto Error;
      }
  
  for (col = 0; col < num_cols; col++)
    {
      for (row = 0; row < num_rows; row++)
        col_vector1[row] = output_block[row][col];
      if (QccVectorInverseDCT(col_vector1,
                              col_vector2,
                              num_rows))
        {
          QccErrorAddMessage("(QccMatrixInverseDCT): Error calling QccVectorInverseDCT()");
          goto Error;
        }
      for (row = 0; row < num_rows; row++)
        output_block[row][col] = col_vector2[row];
    }
  
  return_value = 0;
  goto Return;
 Error:
  return_value = 1;
 Return:
  QccVectorFree(col_vector1);
  QccVectorFree(col_vector2);
  return(return_value);
}


int QccMatrixAddNoiseToRegion(QccMatrix matrix,
                              int num_rows, int num_cols,
                              int start_row, int start_col,
                              int region_num_rows,
                              int region_num_cols,
                              double noise_variance)
{
  int stop_row, stop_col;
  int row, col;
  double val;
  
  if (matrix == NULL)
    return(0);
  
  if ((region_num_rows <= 0) ||
      (region_num_cols <= 0))
    return(0);
  
  if (start_row < 0)
    start_row = 0;
  if (start_col < 0)
    start_col = 0;
  
  stop_row = start_row + region_num_rows - 1;
  stop_col = start_col + region_num_cols - 1;
  
  if (stop_row >= num_rows)
    stop_row = num_rows - 1;
  if (stop_col >= num_cols)
    stop_col = num_cols - 1;
  
  for (row = start_row; row <= stop_row; row++)
    for (col = start_col; col <= stop_col; col++)
      {
        /*
         * The following uses the fact that the variance of
         * a (zero-mean) uniform random variable on the range (-c, c) is
         *     var = (c^2)/3
         */
        
        /*  Uniform(0.0, 1.0), mean = 0.5  */
        val = QccMathRand();
        /*  Uniform(-1.0, 1.0), variance = 1/3, mean = 0.0  */
        val = (val - 0.5)*2.0;
        /*  Uniform(-c, c), Variance = noise_variance  */
        val *= sqrt(3.0*noise_variance);
        
        matrix[row][col] += val;
      }
  
  return(0);
}


int QccMatrixInverse(const QccMatrix matrix,
                     QccMatrix matrix_inverse,
                     int num_rows,
                     int num_cols)
{
#ifdef QCC_USE_GSL

  int return_value;
  gsl_matrix *A = NULL;
  gsl_permutation *P = NULL;
  gsl_matrix *A_inverse = NULL;
  int row, col;
  int signum;
  double determinant;

  if (matrix == NULL)
    return(0);
  if (matrix_inverse == NULL)
    return(0);

  if (num_rows != num_cols)
    {
      QccErrorAddMessage("(QccMatrixInverse): Matrix is not square");
      return(1);
    }

  if ((A = gsl_matrix_alloc(num_rows, num_cols)) == NULL)
    {
      QccErrorAddMessage("(QccMatrixInverse): Error allocating memory");
      goto Error;
    }
  if ((P = gsl_permutation_alloc(num_rows)) == NULL)
    {
      QccErrorAddMessage("(QccMatrixInverse): Error allocating memory");
      goto Error;
    }
  if ((A_inverse = gsl_matrix_alloc(num_rows, num_cols)) == NULL)
    {
      QccErrorAddMessage("(QccMatrixInverse): Error allocating memory");
      goto Error;
    }

  gsl_matrix_set_zero(A);
  for (row = 0; row < num_rows; row++)
    for (col = 0; col < num_cols; col++)
      gsl_matrix_set(A, row, col, matrix[row][col]);

  if (gsl_linalg_LU_decomp(A, P, &signum))
    {
      QccErrorAddMessage("(QccMatrixInverse): Error performing LU decomposition");
      goto Error;
    }

  determinant = gsl_linalg_LU_det(A, signum);
  if (fabs(determinant) < QCCMATHEPS)
    {
      QccErrorAddMessage("(QccMatrixInverse): Matrix must be nonsingular");
      goto Error;
    }

  if (gsl_linalg_LU_invert(A, P, A_inverse))
    {
      QccErrorAddMessage("(QccMatrixInverse): Error performing LU matrix inverse");
      goto Error;
    }
  
  for (row = 0; row < num_rows; row++)
    for (col = 0; col < num_cols; col++)
      matrix_inverse[row][col] = gsl_matrix_get(A_inverse, row, col);
  
  return_value = 0;
  goto Return;
 Error:
  return_value = 1;
 Return:
  gsl_matrix_free(A);
  gsl_permutation_free(P);
  gsl_matrix_free(A_inverse);
  return(return_value);

#else

  QccErrorAddMessage("(QccMatrixInverse): GSL support is not available -- matrix inverse is thus not supported");
  return(1);

#endif
}


int QccMatrixSVD(const QccMatrix A,
                 int num_rows,
                 int num_cols,
                 QccMatrix U,
                 QccVector S,
                 QccMatrix V)
{
#ifdef QCC_USE_GSL

  int return_value;
  gsl_matrix *A2 = NULL;
  gsl_matrix *V2 = NULL;
  gsl_vector *S2 = NULL;
  gsl_vector *work = NULL;
  int row, col;
  int num_rows2;

  if (A == NULL)
    return(0);

  if (num_rows >= num_cols)
    num_rows2 = num_rows;
  else
    num_rows2 = num_cols;

  if ((A2 = gsl_matrix_alloc(num_rows2, num_cols)) == NULL)
    {
      QccErrorAddMessage("(QccMatrixSVD): Error allocating memory");
      goto Error;
    }
  if ((V2 = gsl_matrix_alloc(num_cols, num_cols)) == NULL)
    {
      QccErrorAddMessage("(QccMatrixSVD): Error allocating memory");
      goto Error;
    }
  if ((S2 = gsl_vector_alloc(num_cols)) == NULL)
    {
      QccErrorAddMessage("(QccMatrixSVD): Error allocating memory");
      goto Error;
    }
  if ((work = gsl_vector_alloc(num_cols)) == NULL)
    {
      QccErrorAddMessage("(QccMatrixSVD): Error allocating memory");
      goto Error;
    }
  
  gsl_matrix_set_zero(A2);
  for (row = 0; row < num_rows; row++)
    for (col = 0; col < num_cols; col++)
      gsl_matrix_set(A2, row, col, A[row][col]);
  
  if (gsl_linalg_SV_decomp(A2, V2, S2, work))
    {
      QccErrorAddMessage("(QccMatrixSVD): Error performing SVD");
      goto Error;
    }
  
  if (U != NULL)
    for (row = 0; row < num_rows; row++)
      for (col = 0; col < num_cols; col++)
        U[row][col] = gsl_matrix_get(A2, row, col);
  
  if (S != NULL)
    for (col = 0; col < num_cols; col++)
      S[col] = gsl_vector_get(S2, col);
  
  if (V != NULL)
    for (row = 0; row < num_cols; row++)
      for (col = 0; col < num_cols; col++)
        V[row][col] = gsl_matrix_get(V2, row, col);

  return_value = 0;
  goto Return;
 Error:
  return_value = 1;
 Return:
  gsl_matrix_free(A2);
  gsl_matrix_free(V2);
  gsl_vector_free(S2);
  gsl_vector_free(work);
  return(return_value);

#else

  QccErrorAddMessage("(QccMatrixSVD): GSL support is not available -- SVD is thus not supported");
  return(1);

#endif
}


int QccMatrixOrthogonalize(const QccMatrix matrix1,
                           int num_rows,
                           int num_cols,
                           QccMatrix matrix2,
                           int *num_cols2)
{
  int return_value;
  QccVector eigenvalues = NULL;
  int row, col;
  int num_nonzero_eigenvectors;

  if ((eigenvalues = QccVectorAlloc(num_cols)) == NULL)
    {
      QccErrorAddMessage("(QccMatrixOrthogonalize): Error calling QccVectorAlloc()");
      goto Error;
    }

  if (QccMatrixSVD(matrix1,
                   num_rows,
                   num_cols,
                   matrix2,
                   eigenvalues,
                   NULL))
    {
      QccErrorAddMessage("(QccMatrixOrthogonalize): Error calling QccMatrixSVD()");
      goto Error;
    }

  num_nonzero_eigenvectors = 0;
  for (col = 0; col < num_cols; col++)
    if (eigenvalues[col] < num_rows * eigenvalues[0] * QCCMATHEPS)
      for (row = 0; row < num_rows; row++)
        matrix2[row][col] = 0;
    else
      num_nonzero_eigenvectors++;

  if (num_cols2 != NULL)
    *num_cols2 = num_nonzero_eigenvectors;

  return_value = 0;
  goto Return;
 Error:
  return_value = 1;
 Return:
  QccVectorFree(eigenvalues);
  return(return_value);
}


int QccMatrixNullspace(const QccMatrix matrix1,
                       int num_rows,
                       int num_cols,
                       QccMatrix matrix2,
                       int *num_cols2)
{
  int return_value;
  QccVector eigenvalues = NULL;
  QccMatrix V = NULL;
  int row, col;
  int index;

  if ((eigenvalues = QccVectorAlloc(num_cols)) == NULL)
    {
      QccErrorAddMessage("(QccMatrixOrthogonalize): Error calling QccVectorAlloc()");
      goto Error;
    }
  if ((V = QccMatrixAlloc(num_cols,
                          num_cols)) == NULL)
    {
      QccErrorAddMessage("(QccMatrixOrthogonalize): Error calling QccMatrixAlloc()");
      goto Error;
    }

  if (QccMatrixSVD(matrix1,
                   num_rows,
                   num_cols,
                   NULL,
                   eigenvalues,
                   V))
    {
      QccErrorAddMessage("(QccMatrixOrthogonalize): Error calling QccMatrixSVD()");
      goto Error;
    }

  QccMatrixZero(matrix2, num_cols, num_cols);

  index = 0;
  for (col = 0; col < num_cols; col++)
    if (eigenvalues[col] < num_cols * eigenvalues[0] * QCCMATHEPS)
      {
        for (row = 0; row < num_cols; row++)
          matrix2[row][index] = V[row][col];
        index++;
      }
  
  if (num_cols2 != NULL)
    *num_cols2 = index;

  return_value = 0;
  goto Return;
 Error:
  return_value = 1;
 Return:
  QccVectorFree(eigenvalues);
  QccMatrixFree(V, num_cols);
  return(return_value);
}