Как я могу получить Z*Z^T, используя GSL, где Z - вектор столбца?

Question

Я просматриваю функции GSL для вычисления Z*Z^T, где Z - это n*1 вектор-столбец, но я не смог найти ни одной подходящей функции, любая помощь очень ценится.

Answer 1

GSL поддерживает BLAS (основные подпрограммы линейной алгебры), см. [ http://www.gnu.org/software/gsl/manual/html_node/GSL-BLAS-Interface.html][1].

Функции классифицируются по сложности операции:

• уровень 1: векторно-векторные операции
• уровень 2: матрично-векторные операции
• уровень 3: матрично-матричные операции

Большинство функций имеют разные версии для чисел с плавающей запятой, двойных и комплексных чисел. Ваша операция в основном является внешним произведением вектора Z на себя. Вы можете инициализировать вектор как вектор столбца (здесь числа с двойной точностью):

 gsl_matrix * Z = gsl_matrix_calloc (n,1);

а затем используйте функцию BLAS gsl_blas_dgemm для вычисления Z * Z^T. Первые аргументы этой функции определяют, следует ли транспонировать входные матрицы перед умножением матрицы:

 gsl_blas_dgemm (CblasNoTrans, CblasTrans, 1.0, Z, Z, 0.0, C);

Вот рабочая тестовая программа (вам может понадобиться связать ее с gsl и blas):

#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>

int main(int argc, char ** argv)
{
  size_t n = 4;
  gsl_matrix * Z = gsl_matrix_calloc (n,1);
  gsl_matrix * C = gsl_matrix_calloc (n,n);
  gsl_matrix_set(Z,0,0,1);
  gsl_matrix_set(Z,1,0,2);
  gsl_matrix_set(Z,2,0,0);
  gsl_matrix_set(Z,3,0,1);

  gsl_blas_dgemm (CblasNoTrans,
                  CblasTrans, 1.0, Z, Z, 0.0, C);
  int i,j;
  for (i = 0; i < n; i++)  
    {
      for (j = 0; j < n; j++)
        {
          printf ("%g\t", gsl_matrix_get (C, i, j));
        }
      printf("\n");
    }
  gsl_matrix_free(Z);
  gsl_matrix_free(C);
  return 0;
}