Davood Hajinezhad
Davood Hajinezhad

Reputation: 34

How can I get Z*Z^T using GSL, where Z is column vector?

I am looking through GSL functions to calculate Z*Z^T, where Z is n*1 column vector, but I could not find any fit function, every help is much appreciated.

Upvotes: 1

Views: 594

Answers (1)

StefanW
StefanW

Reputation: 347

GSL supports BLAS (basic linear algebra subprograms), see [http://www.gnu.org/software/gsl/manual/html_node/GSL-BLAS-Interface.html][1].

The functions are classified by the complexity of the operation:

  • level 1: vector-vector operations
  • level 2: matrix-vector operations
  • level 3: matrix-matrix operations

Most functions come in different versions for float, double and complex numbers. Your operation is basically an outer product of the vector Z with itself. You can initialize the vector as a column vector (here double precision numbers):

 gsl_matrix * Z = gsl_matrix_calloc (n,1);

and then use the BLAS function gsl_blas_dgemm to compute Z * Z^T. The first arguments of this function determine, whether or not the input matrices should be transposed before the matrix multiplication:

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

Here's a working test program (you may need to link it against gsl and 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;
}

Upvotes: 1

Related Questions