blas.c

开放gsl矩阵运算

  cblas_ssyrk (CblasRowMajor, Uplo, Trans, INT (N), INT (K), alpha, A->data,
	       INT (A->tda), beta, C->data, INT (C->tda));
  return GSL_SUCCESS;
}


int
gsl_blas_dsyrk (CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t Trans, double alpha,
		const gsl_matrix * A, double beta, gsl_matrix * C)
{
  const size_t M = C->size1;
  const size_t N = C->size2;
  const size_t K = (Trans == CblasNoTrans) ? A->size1 : A->size2;

  if (M != N)
    {
      GSL_ERROR ("matrix C must be square", GSL_ENOTSQR);
    }
  else if (N != K)
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }

  cblas_dsyrk (CblasRowMajor, Uplo, Trans, INT (N), INT (K), alpha, A->data,
	       INT (A->tda), beta, C->data, INT (C->tda));
  return GSL_SUCCESS;

}


int
gsl_blas_csyrk (CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t Trans,
		const gsl_complex_float alpha,
		const gsl_matrix_complex_float * A,
		const gsl_complex_float beta, gsl_matrix_complex_float * C)
{
  const size_t M = C->size1;
  const size_t N = C->size2;
  const size_t MA = (Trans == CblasNoTrans) ? A->size1 : A->size2;

  if (M != N)
    {
      GSL_ERROR ("matrix C must be square", GSL_ENOTSQR);
    }
  else if (N != MA)
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }

  cblas_csyrk (CblasRowMajor, Uplo, Trans, INT (N), INT (MA),
	       GSL_COMPLEX_P (&alpha), A->data, INT (A->tda),
	       GSL_COMPLEX_P (&beta), C->data, INT (C->tda));
  return GSL_SUCCESS;
}


int
gsl_blas_zsyrk (CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t Trans,
		const gsl_complex alpha, const gsl_matrix_complex * A,
		const gsl_complex beta, gsl_matrix_complex * C)
{
  const size_t M = C->size1;
  const size_t N = C->size2;
  const size_t K = (Trans == CblasNoTrans) ? A->size1 : A->size2;

  if (M != N)
    {
      GSL_ERROR ("matrix C must be square", GSL_ENOTSQR);
    }
  else if (N != K)
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }

  cblas_zsyrk (CblasRowMajor, Uplo, Trans, INT (N), INT (K),
	       GSL_COMPLEX_P (&alpha), A->data, INT (A->tda),
	       GSL_COMPLEX_P (&beta), C->data, INT (C->tda));
  return GSL_SUCCESS;
}

/* HERK */

int
gsl_blas_cherk (CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t Trans, float alpha,
		const gsl_matrix_complex_float * A, float beta,
		gsl_matrix_complex_float * C)
{
  const size_t M = C->size1;
  const size_t N = C->size2;
  const size_t K = (Trans == CblasNoTrans) ? A->size1 : A->size2;

  if (M != N)
    {
      GSL_ERROR ("matrix C must be square", GSL_ENOTSQR);
    }
  else if (N != K)
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }

  cblas_cherk (CblasRowMajor, Uplo, Trans, INT (N), INT (K), alpha, A->data,
	       INT (A->tda), beta, C->data, INT (C->tda));
  return GSL_SUCCESS;
}


int
gsl_blas_zherk (CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t Trans, double alpha,
		const gsl_matrix_complex * A, double beta,
		gsl_matrix_complex * C)
{
  const size_t M = C->size1;
  const size_t N = C->size2;
  const size_t K = (Trans == CblasNoTrans) ? A->size1 : A->size2;

  if (M != N)
    {
      GSL_ERROR ("matrix C must be square", GSL_ENOTSQR);
    }
  else if (N != K)
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }

  cblas_zherk (CblasRowMajor, Uplo, Trans, INT (N), INT (K), alpha, A->data,
	       INT (A->tda), beta, C->data, INT (C->tda));
  return GSL_SUCCESS;
}

/* SYR2K */

int
gsl_blas_ssyr2k (CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t Trans, float alpha,
		 const gsl_matrix_float * A, const gsl_matrix_float * B,
		 float beta, gsl_matrix_float * C)
{
  const size_t M = C->size1;
  const size_t N = C->size2;
  const size_t MA = (Trans == CblasNoTrans) ? A->size1 : A->size2;
  const size_t NA = (Trans == CblasNoTrans) ? A->size2 : A->size1;
  const size_t MB = (Trans == CblasNoTrans) ? B->size1 : B->size2;
  const size_t NB = (Trans == CblasNoTrans) ? B->size2 : B->size1;

  if (M != N)
    {
      GSL_ERROR ("matrix C must be square", GSL_ENOTSQR);
    }
  else if (N != MA || N != MB || NA != NB)
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }

  cblas_ssyr2k (CblasRowMajor, Uplo, Trans, INT (N), INT (NA), alpha, A->data,
		INT (A->tda), B->data, INT (B->tda), beta, C->data,
		INT (C->tda));
  return GSL_SUCCESS;
}


int
gsl_blas_dsyr2k (CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t Trans, double alpha,
		 const gsl_matrix * A, const gsl_matrix * B, double beta,
		 gsl_matrix * C)
{
  const size_t M = C->size1;
  const size_t N = C->size2;
  const size_t MA = (Trans == CblasNoTrans) ? A->size1 : A->size2;
  const size_t NA = (Trans == CblasNoTrans) ? A->size2 : A->size1;
  const size_t MB = (Trans == CblasNoTrans) ? B->size1 : B->size2;
  const size_t NB = (Trans == CblasNoTrans) ? B->size2 : B->size1;

  if (M != N)
    {
      GSL_ERROR ("matrix C must be square", GSL_ENOTSQR);
    }
  else if (N != MA || N != MB || NA != NB)
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }

  cblas_dsyr2k (CblasRowMajor, Uplo, Trans, INT (N), INT (NA), alpha, A->data,
		INT (A->tda), B->data, INT (B->tda), beta, C->data,
		INT (C->tda));
  return GSL_SUCCESS;
}


int
gsl_blas_csyr2k (CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t Trans,
		 const gsl_complex_float alpha,
		 const gsl_matrix_complex_float * A,
		 const gsl_matrix_complex_float * B,
		 const gsl_complex_float beta, gsl_matrix_complex_float * C)
{
  const size_t M = C->size1;
  const size_t N = C->size2;
  const size_t MA = (Trans == CblasNoTrans) ? A->size1 : A->size2;
  const size_t NA = (Trans == CblasNoTrans) ? A->size2 : A->size1;
  const size_t MB = (Trans == CblasNoTrans) ? B->size1 : B->size2;
  const size_t NB = (Trans == CblasNoTrans) ? B->size2 : B->size1;

  if (M != N)
    {
      GSL_ERROR ("matrix C must be square", GSL_ENOTSQR);
    }
  else if (N != MA || N != MB || NA != NB)
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }

  cblas_csyr2k (CblasRowMajor, Uplo, Trans, INT (N), INT (NA),
		GSL_COMPLEX_P (&alpha), A->data, INT (A->tda), B->data,
		INT (B->tda), GSL_COMPLEX_P (&beta), C->data, INT (C->tda));
  return GSL_SUCCESS;
}



int
gsl_blas_zsyr2k (CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t Trans,
		 const gsl_complex alpha, const gsl_matrix_complex * A,
		 const gsl_matrix_complex * B, const gsl_complex beta,
		 gsl_matrix_complex * C)
{
  const size_t M = C->size1;
  const size_t N = C->size2;
  const size_t MA = (Trans == CblasNoTrans) ? A->size1 : A->size2;
  const size_t NA = (Trans == CblasNoTrans) ? A->size2 : A->size1;
  const size_t MB = (Trans == CblasNoTrans) ? B->size1 : B->size2;
  const size_t NB = (Trans == CblasNoTrans) ? B->size2 : B->size1;

  if (M != N)
    {
      GSL_ERROR ("matrix C must be square", GSL_ENOTSQR);
    }
  else if (N != MA || N != MB || NA != NB)
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }

  cblas_zsyr2k (CblasRowMajor, Uplo, Trans, INT (N), INT (NA),
		GSL_COMPLEX_P (&alpha), A->data, INT (A->tda), B->data,
		INT (B->tda), GSL_COMPLEX_P (&beta), C->data, INT (C->tda));
  return GSL_SUCCESS;
}

/* HER2K */

int
gsl_blas_cher2k (CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t Trans,
		 const gsl_complex_float alpha,
		 const gsl_matrix_complex_float * A,
		 const gsl_matrix_complex_float * B, float beta,
		 gsl_matrix_complex_float * C)
{
  const size_t M = C->size1;
  const size_t N = C->size2;
  const size_t MA = (Trans == CblasNoTrans) ? A->size1 : A->size2;
  const size_t NA = (Trans == CblasNoTrans) ? A->size2 : A->size1;
  const size_t MB = (Trans == CblasNoTrans) ? B->size1 : B->size2;
  const size_t NB = (Trans == CblasNoTrans) ? B->size2 : B->size1;

  if (M != N)
    {
      GSL_ERROR ("matrix C must be square", GSL_ENOTSQR);
    }
  else if (N != MA || N != MB || NA != NB)
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }

  cblas_cher2k (CblasRowMajor, Uplo, Trans, INT (N), INT (NA),
		GSL_COMPLEX_P (&alpha), A->data, INT (A->tda), B->data,
		INT (B->tda), beta, C->data, INT (C->tda));
  return GSL_SUCCESS;

}


int
gsl_blas_zher2k (CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t Trans,
		 const gsl_complex alpha, const gsl_matrix_complex * A,
		 const gsl_matrix_complex * B, double beta,
		 gsl_matrix_complex * C)
{
  const size_t M = C->size1;
  const size_t N = C->size2;
  const size_t MA = (Trans == CblasNoTrans) ? A->size1 : A->size2;
  const size_t NA = (Trans == CblasNoTrans) ? A->size2 : A->size1;
  const size_t MB = (Trans == CblasNoTrans) ? B->size1 : B->size2;
  const size_t NB = (Trans == CblasNoTrans) ? B->size2 : B->size1;

  if (M != N)
    {
      GSL_ERROR ("matrix C must be square", GSL_ENOTSQR);
    }
  else if (N != MA || N != MB || NA != NB)
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }

  cblas_zher2k (CblasRowMajor, Uplo, Trans, INT (N), INT (NA),
		GSL_COMPLEX_P (&alpha), A->data, INT (A->tda), B->data,
		INT (B->tda), beta, C->data, INT (C->tda));
  return GSL_SUCCESS;

}

/* TRMM */

int
gsl_blas_strmm (CBLAS_SIDE_t Side, CBLAS_UPLO_t Uplo,
		CBLAS_TRANSPOSE_t TransA, CBLAS_DIAG_t Diag, float alpha,
		const gsl_matrix_float * A, gsl_matrix_float * B)
{
  const size_t M = B->size1;
  const size_t N = B->size2;
  const size_t MA = A->size1;
  const size_t NA = A->size2;

  if (MA != NA)
    {
      GSL_ERROR ("matrix A must be square", GSL_ENOTSQR);
    }

  if ((Side == CblasLeft && M == MA) || (Side == CblasRight && N == MA))
    {
      cblas_strmm (CblasRowMajor, Side, Uplo, TransA, Diag, INT (M), INT (N),
		   alpha, A->data, INT (A->tda), B->data, INT (B->tda));
      return GSL_SUCCESS;
    }
  else
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }
}


int
gsl_blas_dtrmm (CBLAS_SIDE_t Side, CBLAS_UPLO_t Uplo,
		CBLAS_TRANSPOSE_t TransA, CBLAS_DIAG_t Diag, double alpha,
		const gsl_matrix * A, gsl_matrix * B)
{
  const size_t M = B->size1;
  const size_t N = B->size2;
  const size_t MA = A->size1;
  const size_t NA = A->size2;

  if (MA != NA)
    {
      GSL_ERROR ("matrix A must be square", GSL_ENOTSQR);
    }

  if ((Side == CblasLeft && M == MA) || (Side == CblasRight && N == MA))
    {
      cblas_dtrmm (CblasRowMajor, Side, Uplo, TransA, Diag, INT (M), INT (N),
		   alpha, A->data, INT (A->tda), B->data, INT (B->tda));
      return GSL_SUCCESS;
    }
  else
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }
}


int
gsl_blas_ctrmm (CBLAS_SIDE_t Side, CBLAS_UPLO_t Uplo,
		CBLAS_TRANSPOSE_t TransA, CBLAS_DIAG_t Diag,
		const gsl_complex_float alpha,
		const gsl_matrix_complex_float * A,
		gsl_matrix_complex_float * B)
{
  const size_t M = B->size1;
  const size_t N = B->size2;
  const size_t MA = A->size1;
  const size_t NA = A->size2;

  if (MA != NA)
    {
      GSL_ERROR ("matrix A must be square", GSL_ENOTSQR);
    }

  if ((Side == CblasLeft && M == MA) || (Side == CblasRight && N == MA))
    {
      cblas_ctrmm (CblasRowMajor, Side, Uplo, TransA, Diag, INT (M), INT (N),
		   GSL_COMPLEX_P (&alpha), A->data, INT (A->tda), B->data,
		   INT (B->tda));
      return GSL_SUCCESS;
    }
  else
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }
}


int
gsl_blas_ztrmm (CBLAS_SIDE_t Side, CBLAS_UPLO_t Uplo,
		CBLAS_TRANSPOSE_t TransA, CBLAS_DIAG_t Diag,
		const gsl_complex alpha, const gsl_matrix_complex * A,
		gsl_matrix_complex * B)
{
  const size_t M = B->size1;
  const size_t N = B->size2;
  const size_t MA = A->size1;
  const size_t NA = A->size2;

  if (MA != NA)
    {
      GSL_ERROR ("matrix A must be square", GSL_ENOTSQR);
    }

  if ((Side == CblasLeft && M == MA) || (Side == CblasRight && N == MA))
    {
      cblas_ztrmm (CblasRowMajor, Side, Uplo, TransA, Diag, INT (M), INT (N),
		   GSL_COMPLEX_P (&alpha), A->data, INT (A->tda), B->data,
		   INT (B->tda));
      return GSL_SUCCESS;
    }
  else
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }
}


/* TRSM */

int
gsl_blas_strsm (CBLAS_SIDE_t Side, CBLAS_UPLO_t Uplo,
		CBLAS_TRANSPOSE_t TransA, CBLAS_DIAG_t Diag, float alpha,
		const gsl_matrix_float * A, gsl_matrix_float * B)
{
  const size_t M = B->size1;
  const size_t N = B->size2;
  const size_t MA = A->size1;
  const size_t NA = A->size2;

  if (MA != NA)
    {
      GSL_ERROR ("matrix A must be square", GSL_ENOTSQR);
    }

  if ((Side == CblasLeft && M == MA) || (Side == CblasRight && N == MA))
    {
      cblas_strsm (CblasRowMajor, Side, Uplo, TransA, Diag, INT (M), INT (N),
		   alpha, A->data, INT (A->tda), B->data, INT (B->tda));
      return GSL_SUCCESS;
    }
  else
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }
}


int
gsl_blas_dtrsm (CBLAS_SIDE_t Side, CBLAS_UPLO_t Uplo,
		CBLAS_TRANSPOSE_t TransA, CBLAS_DIAG_t Diag, double alpha,
		const gsl_matrix * A, gsl_matrix * B)
{
  const size_t M = B->size1;
  const size_t N = B->size2;
  const size_t MA = A->size1;
  const size_t NA = A->size2;

  if (MA != NA)
    {
      GSL_ERROR ("matrix A must be square", GSL_ENOTSQR);
    }

  if ((Side == CblasLeft && M == MA) || (Side == CblasRight && N == MA))
    {
      cblas_dtrsm (CblasRowMajor, Side, Uplo, TransA, Diag, INT (M), INT (N),
		   alpha, A->data, INT (A->tda), B->data, INT (B->tda));
      return GSL_SUCCESS;
    }
  else
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }
}


int
gsl_blas_ctrsm (CBLAS_SIDE_t Side, CBLAS_UPLO_t Uplo,
		CBLAS_TRANSPOSE_t TransA, CBLAS_DIAG_t Diag,
		const gsl_complex_float alpha,
		const gsl_matrix_complex_float * A,
		gsl_matrix_complex_float * B)
{
  const size_t M = B->size1;
  const size_t N = B->size2;
  const size_t MA = A->size1;
  const size_t NA = A->size2;

  if (MA != NA)
    {
      GSL_ERROR ("matrix A must be square", GSL_ENOTSQR);
    }

  if ((Side == CblasLeft && M == MA) || (Side == CblasRight && N == MA))
    {
      cblas_ctrsm (CblasRowMajor, Side, Uplo, TransA, Diag, INT (M), INT (N),
		   GSL_COMPLEX_P (&alpha), A->data, INT (A->tda), B->data,
		   INT (B->tda));
      return GSL_SUCCESS;
    }
  else
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }
}


int
gsl_blas_ztrsm (CBLAS_SIDE_t Side, CBLAS_UPLO_t Uplo,
		CBLAS_TRANSPOSE_t TransA, CBLAS_DIAG_t Diag,
		const gsl_complex alpha, const gsl_matrix_complex * A,
		gsl_matrix_complex * B)
{
  const size_t M = B->size1;
  const size_t N = B->size2;
  const size_t MA = A->size1;
  const size_t NA = A->size2;

  if (MA != NA)
    {
      GSL_ERROR ("matrix A must be square", GSL_ENOTSQR);
    }

  if ((Side == CblasLeft && M == MA) || (Side == CblasRight && N == MA))
    {
      cblas_ztrsm (CblasRowMajor, Side, Uplo, TransA, Diag, INT (M), INT (N),
		   GSL_COMPLEX_P (&alpha), A->data, INT (A->tda), B->data,
		   INT (B->tda));
      return GSL_SUCCESS;
    }
  else
    {
      GSL_ERROR ("invalid length", GSL_EBADLEN);
    }
}