Search code examples
c++gsl

best way to get outer product in GSL

I have two vectors, x, y and I need to get outer product i.e xy'. I have to do this in GSL library and I am not sure if there is a function to do it.

Can anyone tell me if there is a function or procedure to calculate the outer product in GSL?

Solution

  • I believe the function you are looking for is in the BLAS interface, which has the following prototype:

    int gsl_blas_dsyr(CBLAS_UPLO_t Uplo, double alpha, const gsl_vector * x, gsl_matrix * A

    This function sends A -> A + alpha x x'. To calculate the outer product, take the case alpha = 1, and A=0. Note, that Uplo = CblasUpper specifies that, since A is symmetric, only the upper-triangular entries are computed. If you want the full matrix, be sure to symmeterize after. So,

    int outerProduct(gsl_matrix* A, const gsl_vector *x, const int len){
  A = gsl_matrix_calloc(len,len);
  return gsl_blas_dsyr(CblasUpper,1.0,x,A);
}

    will allocate a len by len matrix, and, assuming x is of length len, will produce the outer-product into the upper-triangular portion of A.