[Top][All Lists]
[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;
}
- [Help-gsl] undefined reference to gsl_matrix_complex_const_subrow,
Jo_Jo <=