help-gsl
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

[Help-gsl] undefined reference to gsl_matrix_complex_const_subrow


From: Jo_Jo
Subject: [Help-gsl] undefined reference to gsl_matrix_complex_const_subrow
Date: Wed, 05 Dec 2007 17:46:42 +0100
User-agent: Thunderbird 2.0.0.6 (X11/20070801)

Hi gsl-help,

In my compiler line command:
"gcc -Wall cholesky_complex.c -lgsl -lgslcblas -lm -o cholesky_complex"
become I in output:
"
/tmp/ccntKlUB.o: In function `gsl_linalg_complex_cholesky_decomp':
cholesky_complex.c:(.text+0xd6): undefined reference to
`gsl_matrix_complex_const_subrow'
cholesky_complex.c:(.text+0x1f3): undefined reference to
`gsl_matrix_complex_subcolumn'
cholesky_complex.c:(.text+0x22b): undefined reference to
`gsl_matrix_complex_subrow'
collect2: ld returned 1 exit status "

Where is mistake? Please run my program


// program cholesky_complex.c
#include <stdio.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_permutation.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_matrix_complex_double.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_errno.h>

/*--------------------------------------------------------*/

/* linalg/choleskyc.c
 *
 * Copyright (C) 2007 Patrick Alken
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
02110-1301, USA.
 */
/*
#include <config.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_errno.h> */

/*
 * This module contains routines related to the Cholesky decomposition
 * of a complex Hermitian positive definite matrix.
 */

static void cholesky_complex_conj_vector(gsl_vector_complex *v);

/*
gsl_linalg_complex_cholesky_decomp()
  Perform the Cholesky decomposition on a Hermitian positive definite
matrix. See Golub & Van Loan, "Matrix Computations" (3rd ed),
algorithm 4.2.2.

Inputs: A - (input/output) complex postive definite matrix

Return: success or error

The lower triangle of A is overwritten with the Cholesky decomposition
*/

int
gsl_linalg_complex_cholesky_decomp(gsl_matrix_complex *A)
{
  const size_t N = A->size1;
 
  if (N != A->size2)
    {
      GSL_ERROR("cholesky decomposition requires square matrix",
GSL_ENOTSQR);
    }
  else
    {
      size_t j;
      gsl_complex z;
      double ajj;

      for (j = 0; j < N; ++j)
        {
          z = gsl_matrix_complex_get(A, j, j);
          ajj = GSL_REAL(z);

          if (j > 0)
            {
              gsl_vector_complex_const_view aj =
                gsl_matrix_complex_const_subrow(A, j, 0, j);

              gsl_blas_zdotc(&aj.vector, &aj.vector, &z);
              ajj -= GSL_REAL(z);
            }

          if (ajj <= 0.0)
            {
              GSL_ERROR("matrix is not positive definite", GSL_EDOM);
            }

          ajj = sqrt(ajj);
          GSL_SET_COMPLEX(&z, ajj, 0.0);
          gsl_matrix_complex_set(A, j, j, z);

          if (j < N - 1)
            {
              gsl_vector_complex_view av =
                gsl_matrix_complex_subcolumn(A, j, j + 1, N - j - 1);

              if (j > 0)
                {
                  gsl_vector_complex_view aj =
                    gsl_matrix_complex_subrow(A, j, 0, j);
                  gsl_matrix_complex_view am =
                    gsl_matrix_complex_submatrix(A, j + 1, 0, N - j - 1, j);

                  cholesky_complex_conj_vector(&aj.vector);

                  gsl_blas_zgemv(CblasNoTrans,
                                 GSL_COMPLEX_NEGONE,
                                 &am.matrix,
                                 &aj.vector,
                                 GSL_COMPLEX_ONE,
                                 &av.vector);

                  cholesky_complex_conj_vector(&aj.vector);
                }

              gsl_blas_zdscal(1.0 / ajj, &av.vector);
            }
        }

      return GSL_SUCCESS;
    }
} /* gsl_linalg_complex_cholesky_decomp() */

/*
gsl_linalg_complex_cholesky_solve()
  Solve A x = b where A is in cholesky form
*/

int
gsl_linalg_complex_cholesky_solve (const gsl_matrix_complex * cholesky,
                                   const gsl_vector_complex * b,
                                   gsl_vector_complex * x)
{
  if (cholesky->size1 != cholesky->size2)
    {
      GSL_ERROR ("cholesky matrix must be square", GSL_ENOTSQR);
    }
  else if (cholesky->size1 != b->size)
    {
      GSL_ERROR ("matrix size must match b size", GSL_EBADLEN);
    }
  else if (cholesky->size2 != x->size)
    {
      GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
    }
  else
    {
      gsl_vector_complex_memcpy (x, b);

      /* solve for y using forward-substitution, L y = b */

      gsl_blas_ztrsv (CblasLower, CblasNoTrans, CblasNonUnit, cholesky, x);

      /* perform back-substitution, L^H x = y */

      gsl_blas_ztrsv (CblasLower, CblasConjTrans, CblasNonUnit,
cholesky, x);

      return GSL_SUCCESS;
    }
} /* gsl_linalg_complex_cholesky_solve() */

/*
gsl_linalg_complex_cholesky_svx()
  Solve A x = b in place where A is in cholesky form
*/

int
gsl_linalg_complex_cholesky_svx (const gsl_matrix_complex * cholesky,
                                 gsl_vector_complex * x)
{
  if (cholesky->size1 != cholesky->size2)
    {
      GSL_ERROR ("cholesky matrix must be square", GSL_ENOTSQR);
    }
  else if (cholesky->size2 != x->size)
    {
      GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
    }
  else
    {
      /* solve for y using forward-substitution, L y = b */

      gsl_blas_ztrsv (CblasLower, CblasNoTrans, CblasNonUnit, cholesky, x);

      /* perform back-substitution, L^H x = y */

      gsl_blas_ztrsv (CblasLower, CblasConjTrans, CblasNonUnit,
cholesky, x);

      return GSL_SUCCESS;
    }
} /* gsl_linalg_complex_cholesky_svx() */

/********************************************
 *           INTERNAL ROUTINES              *
 ********************************************/

static void
cholesky_complex_conj_vector(gsl_vector_complex *v)
{
  size_t i;

  for (i = 0; i < v->size; ++i)
    {
      gsl_complex z = gsl_vector_complex_get(v, i);
      gsl_vector_complex_set(v, i, gsl_complex_conjugate(z));
    }
} /* cholesky_complex_conj_vector() */


/* end of linalg/choleskyc.c */
/*-----------------------------------------------------------*/

int main(){ int M = 4, N = 4, i, j;
double re_data[] = {  4.18,  0.60,  0.57,  0.96,
                      0.60,  5.28,  0.30,  0.58,
                      0.57,  0.30,  6.97,  0.19,
                      0.96,  0.58,  0.19,  7.85 };
double imag_data[] = {  0.18,  0.69,  0.57,  0.60,
                        0.69,  0.28,  0.91,  0.80,
                        0.57,  0.91,  0.97,  0.19,
                        0.60,  0.80,  0.19,  0.85 };
double bz_real_data[] = { 1.1, 2.2, 3.3, 4.4 };
double bz_imag_data[]  = { 5.5, 6.6, 7.7, 8.8 };
printf("\n  copy matrix_view to matrix ----------------\n ");
double temp1;
gsl_matrix_view Re_v = gsl_matrix_view_array (re_data, M, N);
gsl_matrix *Re = gsl_matrix_alloc ( M, N );
// copy Re_v -> Re  
   for (i = 0; i < M; i++)
       for (j = 0; j < N; j++){
           temp1 = gsl_matrix_get (&Re_v.matrix, i, j);
        gsl_matrix_set(Re, i, j, temp1) ; }   
gsl_matrix_view Im_v = gsl_matrix_view_array (imag_data, M, N);
gsl_matrix *Im = gsl_matrix_alloc ( M, N );
// copy Im_v -> Im  
   for (i = 0; i < M; i++)
       for (j = 0; j < N; j++){
           temp1 = gsl_matrix_get (&Im_v.matrix, i, j);
        gsl_matrix_set(Im, i, j, temp1) ; }   
// complex matrix Z       
gsl_matrix_complex *Z= gsl_matrix_complex_calloc (M, N);
gsl_complex temp;
double temp2;
for (i = 0; i < M; i++)
       for (j = 0; j < N; j++){
           temp1 = gsl_matrix_get (Re, i, j);
           temp2 = gsl_matrix_get (Im, i, j);
           temp = gsl_complex_rect (temp1, temp2 );
        gsl_matrix_complex_set(Z, i, j, temp) ;};
// complex vector b_z
 gsl_vector_complex * b_z  = gsl_vector_complex_alloc(M);
  for(i=0;i<M;i++){
       temp1= bz_real_data[i];
       temp2= bz_imag_data[i];
       temp = gsl_complex_rect(temp1, temp2);
       gsl_vector_complex_set(b_z,i,temp); };
//complex vector x_z
  gsl_vector_complex *x_z= gsl_vector_complex_alloc(M);
//vector x_z
(void) gsl_linalg_complex_cholesky_decomp ( Z ) ;
(void)  gsl_linalg_complex_cholesky_solve ( Z, b_z, x_z ) ;
(void)gsl_vector_complex_fprintf (stdout, x_z, "%f   %f" );

gsl_matrix_free(Re);
gsl_matrix_free(Im);
gsl_matrix_complex_free(Z);

 gsl_vector_complex_free (x_z);
 gsl_vector_complex_free (b_z);
      
 return 0;
}
 


Best regards

Andrzej Buczynski
Moosburg, Germany



// program cholesky_complex.c
#include <stdio.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_permutation.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_matrix_complex_double.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_errno.h>

/*--------------------------------------------------------*/

/* linalg/choleskyc.c
 * 
 * Copyright (C) 2007 Patrick Alken
 * 
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3 of the License, or (at
 * your option) any later version.
 * 
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, 
USA.
 */
/*
#include <config.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_errno.h> */

/*
 * This module contains routines related to the Cholesky decomposition
 * of a complex Hermitian positive definite matrix.
 */

static void cholesky_complex_conj_vector(gsl_vector_complex *v);

/*
gsl_linalg_complex_cholesky_decomp()
  Perform the Cholesky decomposition on a Hermitian positive definite
matrix. See Golub & Van Loan, "Matrix Computations" (3rd ed),
algorithm 4.2.2.

Inputs: A - (input/output) complex postive definite matrix

Return: success or error

The lower triangle of A is overwritten with the Cholesky decomposition
*/

int
gsl_linalg_complex_cholesky_decomp(gsl_matrix_complex *A)
{
  const size_t N = A->size1;
  
  if (N != A->size2)
    {
      GSL_ERROR("cholesky decomposition requires square matrix", GSL_ENOTSQR);
    }
  else
    {
      size_t j;
      gsl_complex z;
      double ajj;

      for (j = 0; j < N; ++j)
        {
          z = gsl_matrix_complex_get(A, j, j);
          ajj = GSL_REAL(z);

          if (j > 0)
            {
              gsl_vector_complex_const_view aj =
                gsl_matrix_complex_const_subrow(A, j, 0, j);

              gsl_blas_zdotc(&aj.vector, &aj.vector, &z);
              ajj -= GSL_REAL(z);
            }

          if (ajj <= 0.0)
            {
              GSL_ERROR("matrix is not positive definite", GSL_EDOM);
            }

          ajj = sqrt(ajj);
          GSL_SET_COMPLEX(&z, ajj, 0.0);
          gsl_matrix_complex_set(A, j, j, z);

          if (j < N - 1)
            {
              gsl_vector_complex_view av =
                gsl_matrix_complex_subcolumn(A, j, j + 1, N - j - 1);

              if (j > 0)
                {
                  gsl_vector_complex_view aj =
                    gsl_matrix_complex_subrow(A, j, 0, j);
                  gsl_matrix_complex_view am =
                    gsl_matrix_complex_submatrix(A, j + 1, 0, N - j - 1, j);

                  cholesky_complex_conj_vector(&aj.vector);

                  gsl_blas_zgemv(CblasNoTrans,
                                 GSL_COMPLEX_NEGONE,
                                 &am.matrix,
                                 &aj.vector,
                                 GSL_COMPLEX_ONE,
                                 &av.vector);

                  cholesky_complex_conj_vector(&aj.vector);
                }

              gsl_blas_zdscal(1.0 / ajj, &av.vector);
            }
        }

      return GSL_SUCCESS;
    }
} /* gsl_linalg_complex_cholesky_decomp() */

/*
gsl_linalg_complex_cholesky_solve()
  Solve A x = b where A is in cholesky form
*/

int
gsl_linalg_complex_cholesky_solve (const gsl_matrix_complex * cholesky,
                                   const gsl_vector_complex * b,
                                   gsl_vector_complex * x)
{
  if (cholesky->size1 != cholesky->size2)
    {
      GSL_ERROR ("cholesky matrix must be square", GSL_ENOTSQR);
    }
  else if (cholesky->size1 != b->size)
    {
      GSL_ERROR ("matrix size must match b size", GSL_EBADLEN);
    }
  else if (cholesky->size2 != x->size)
    {
      GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
    }
  else
    {
      gsl_vector_complex_memcpy (x, b);

      /* solve for y using forward-substitution, L y = b */

      gsl_blas_ztrsv (CblasLower, CblasNoTrans, CblasNonUnit, cholesky, x);

      /* perform back-substitution, L^H x = y */

      gsl_blas_ztrsv (CblasLower, CblasConjTrans, CblasNonUnit, cholesky, x);

      return GSL_SUCCESS;
    }
} /* gsl_linalg_complex_cholesky_solve() */

/*
gsl_linalg_complex_cholesky_svx()
  Solve A x = b in place where A is in cholesky form
*/

int
gsl_linalg_complex_cholesky_svx (const gsl_matrix_complex * cholesky,
                                 gsl_vector_complex * x)
{
  if (cholesky->size1 != cholesky->size2)
    {
      GSL_ERROR ("cholesky matrix must be square", GSL_ENOTSQR);
    }
  else if (cholesky->size2 != x->size)
    {
      GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
    }
  else
    {
      /* solve for y using forward-substitution, L y = b */

      gsl_blas_ztrsv (CblasLower, CblasNoTrans, CblasNonUnit, cholesky, x);

      /* perform back-substitution, L^H x = y */

      gsl_blas_ztrsv (CblasLower, CblasConjTrans, CblasNonUnit, cholesky, x);

      return GSL_SUCCESS;
    }
} /* gsl_linalg_complex_cholesky_svx() */

/********************************************
 *           INTERNAL ROUTINES              *
 ********************************************/

static void
cholesky_complex_conj_vector(gsl_vector_complex *v)
{
  size_t i;

  for (i = 0; i < v->size; ++i)
    {
      gsl_complex z = gsl_vector_complex_get(v, i);
      gsl_vector_complex_set(v, i, gsl_complex_conjugate(z));
    }
} /* cholesky_complex_conj_vector() */


/* end of linalg/choleskyc.c */
/*-----------------------------------------------------------*/

int main(){ int M = 4, N = 4, i, j;
double re_data[] = {  4.18,  0.60,  0.57,  0.96,
                      0.60,  5.28,  0.30,  0.58,
                      0.57,  0.30,  6.97,  0.19,
                      0.96,  0.58,  0.19,  7.85 }; 
double imag_data[] = {  0.18,  0.69,  0.57,  0.60,
                        0.69,  0.28,  0.91,  0.80,
                        0.57,  0.91,  0.97,  0.19,
                        0.60,  0.80,  0.19,  0.85 }; 
double bz_real_data[] = { 1.1, 2.2, 3.3, 4.4 };
double bz_imag_data[]  = { 5.5, 6.6, 7.7, 8.8 };
printf("\n  copy matrix_view to matrix ----------------\n ");
double temp1;
gsl_matrix_view Re_v = gsl_matrix_view_array (re_data, M, N);
gsl_matrix *Re = gsl_matrix_alloc ( M, N );
// copy Re_v -> Re   
   for (i = 0; i < M; i++)
       for (j = 0; j < N; j++){
        temp1 = gsl_matrix_get (&Re_v.matrix, i, j);
        gsl_matrix_set(Re, i, j, temp1) ; }     
gsl_matrix_view Im_v = gsl_matrix_view_array (imag_data, M, N);
gsl_matrix *Im = gsl_matrix_alloc ( M, N );
// copy Im_v -> Im   
   for (i = 0; i < M; i++)
       for (j = 0; j < N; j++){
        temp1 = gsl_matrix_get (&Im_v.matrix, i, j);
        gsl_matrix_set(Im, i, j, temp1) ; }     
// complex matrix Z        
gsl_matrix_complex *Z= gsl_matrix_complex_calloc (M, N);
gsl_complex temp;
double temp2;
for (i = 0; i < M; i++)
       for (j = 0; j < N; j++){
        temp1 = gsl_matrix_get (Re, i, j);
        temp2 = gsl_matrix_get (Im, i, j);
        temp = gsl_complex_rect (temp1, temp2 );
        gsl_matrix_complex_set(Z, i, j, temp) ;};
// complex vector b_z
 gsl_vector_complex * b_z  = gsl_vector_complex_alloc(M);
  for(i=0;i<M;i++){
         temp1= bz_real_data[i];
         temp2= bz_imag_data[i]; 
         temp = gsl_complex_rect(temp1, temp2); 
         gsl_vector_complex_set(b_z,i,temp); };
//complex vector x_z
  gsl_vector_complex *x_z= gsl_vector_complex_alloc(M);
//vector x_z
(void) gsl_linalg_complex_cholesky_decomp ( Z ) ;
(void)  gsl_linalg_complex_cholesky_solve ( Z, b_z, x_z ) ;
(void)gsl_vector_complex_fprintf (stdout, x_z, "%f   %f" );

gsl_matrix_free(Re);
gsl_matrix_free(Im);
gsl_matrix_complex_free(Z);

 gsl_vector_complex_free (x_z);
 gsl_vector_complex_free (b_z);
       
 return 0;
}
 

reply via email to

[Prev in Thread] Current Thread [Next in Thread]