/*
blas2.c
Ruby/GSL: Ruby extension library for GSL (GNU Scientific Library)
(C) Copyright 2001-2006 by Yoshiki Tsunesada
Ruby/GSL is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License.
This library is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY.
*/
#include "rb_gsl_config.h"
#include <gsl/gsl_blas.h>
#include "rb_gsl_common.h"
#include "rb_gsl_complex.h"
#include "rb_gsl_array.h"
static void define_const(VALUE module);
static void define_const(VALUE module)
{
rb_define_const(module, "CblasRowMajor", INT2FIX(CblasRowMajor));
rb_define_const(module, "CblasColMajor", INT2FIX(CblasColMajor));
rb_define_const(module, "RowMajor", INT2FIX(CblasRowMajor));
rb_define_const(module, "ColMajor", INT2FIX(CblasColMajor));
rb_define_const(module, "CblasNoTrans", INT2FIX(CblasNoTrans));
rb_define_const(module, "CblasTrans", INT2FIX(CblasTrans));
rb_define_const(module, "CblasConjTrans", INT2FIX(CblasConjTrans));
rb_define_const(module, "NoTrans", INT2FIX(CblasNoTrans));
rb_define_const(module, "Trans", INT2FIX(CblasTrans));
rb_define_const(module, "ConjTrans", INT2FIX(CblasConjTrans));
rb_define_const(module, "CblasUpper", INT2FIX(CblasUpper));
rb_define_const(module, "CblasLower", INT2FIX(CblasLower));
rb_define_const(module, "Upper", INT2FIX(CblasUpper));
rb_define_const(module, "Lower", INT2FIX(CblasLower));
rb_define_const(module, "CblasNonUnit", INT2FIX(CblasNonUnit));
rb_define_const(module, "CblasUnit", INT2FIX(CblasUnit));
rb_define_const(module, "NonUnit", INT2FIX(CblasNonUnit));
rb_define_const(module, "Unit", INT2FIX(CblasUnit));
rb_define_const(module, "CblasLeft", INT2FIX(CblasLeft));
rb_define_const(module, "CblasRight", INT2FIX(CblasRight));
rb_define_const(module, "Left", INT2FIX(CblasLeft));
rb_define_const(module, "Right", INT2FIX(CblasRight));
}
static VALUE rb_gsl_blas_dgemv(int argc, VALUE *argv, VALUE obj)
{
gsl_matrix *A = NULL;
gsl_vector *x = NULL, *y = NULL;
double a, b;
int type, istart, flag = 0;
switch (TYPE(obj)) {
case T_MODULE:
case T_CLASS:
case T_OBJECT:
if (argc < 4) rb_raise(rb_eArgError, "wrong number of arguments (%d for >= 4)",
argc);
CHECK_FIXNUM(argv[0]);
Need_Float(argv[1]);
CHECK_MATRIX(argv[2]);
CHECK_VECTOR(argv[3]);
type = FIX2INT(argv[0]);
a = NUM2DBL(argv[1]);
Data_Get_Struct(argv[2], gsl_matrix, A);
Data_Get_Struct(argv[3], gsl_vector, x);
istart = 4;
break;
default:
Data_Get_Struct(obj, gsl_matrix, A);
if (argc < 3) rb_raise(rb_eArgError, "wrong number of arguments (%d for >= 3)",
argc);
CHECK_FIXNUM(argv[0]);
Need_Float(argv[1]);
CHECK_VECTOR(argv[2]);
type = FIX2INT(argv[0]);
a = NUM2DBL(argv[1]);
Data_Get_Struct(argv[2], gsl_vector, x);
istart = 3;
break;
}
switch (argc - istart) {
case 2:
Need_Float(argv[istart]);
CHECK_VECTOR(argv[istart+1]);
b = NUM2DBL(argv[istart]);
Data_Get_Struct(argv[istart+1], gsl_vector, y);
break;
case 0:
b = 0.0;
y = gsl_vector_alloc(x->size);
flag = 1;
break;
default:
rb_raise(rb_eArgError, "wrong number of arguments");
break;
}
gsl_blas_dgemv(type, a, A, x, b, y);
if (flag == 1) return Data_Wrap_Struct(cgsl_vector, 0, gsl_vector_free, y);
else return argv[argc-1];
}
static VALUE rb_gsl_blas_dgemv2(int argc, VALUE *argv, VALUE obj)
{
gsl_matrix *A = NULL;
gsl_vector *x = NULL, *y, *ynew;
double a, b;
int type, istart, flag = 0;
switch (TYPE(obj)) {
case T_MODULE:
case T_CLASS:
case T_OBJECT:
if (argc < 4) rb_raise(rb_eArgError, "wrong number of arguments (%d for >= 4)",
argc);
CHECK_FIXNUM(argv[0]);
Need_Float(argv[1]);
CHECK_MATRIX(argv[2]);
CHECK_VECTOR(argv[3]);
type = FIX2INT(argv[0]);
a = NUM2DBL(argv[1]);
Data_Get_Struct(argv[2], gsl_matrix, A);
Data_Get_Struct(argv[3], gsl_vector, x);
istart = 4;
break;
default:
Data_Get_Struct(obj, gsl_matrix, A);
if (argc < 3) rb_raise(rb_eArgError, "wrong number of arguments (%d for >= 3)",
argc);
CHECK_FIXNUM(argv[0]);
Need_Float(argv[1]);
CHECK_VECTOR(argv[2]);
type = FIX2INT(argv[0]);
a = NUM2DBL(argv[1]);
Data_Get_Struct(argv[2], gsl_vector, x);
istart = 3;
break;
}
switch (argc - istart) {
case 2:
Need_Float(argv[istart]);
CHECK_VECTOR(argv[istart+1]);
b = NUM2DBL(argv[istart]);
Data_Get_Struct(argv[istart+1], gsl_vector, y);
break;
case 0:
b = 0.0;
y = gsl_vector_alloc(x->size);
flag = 1;
break;
default:
rb_raise(rb_eArgError, "wrong number of arguments");
break;
}
ynew = gsl_vector_alloc(y->size);
gsl_vector_memcpy(ynew, y);
gsl_blas_dgemv(type, a, A, x, b, ynew);
if (flag == 1) gsl_vector_free(y);
return Data_Wrap_Struct(cgsl_vector, 0, gsl_vector_free, ynew);
}
static VALUE rb_gsl_blas_zgemv(int argc, VALUE *argv, VALUE obj)
{
gsl_matrix_complex *A = NULL;
gsl_vector_complex *x = NULL, *y = NULL;
gsl_complex *a, *b, z;
int type, istart, flag = 0;
switch (TYPE(obj)) {
case T_MODULE:
case T_CLASS:
case T_OBJECT:
if (argc < 4) rb_raise(rb_eArgError, "wrong number of arguments (%d for >= 4)",
argc);
CHECK_FIXNUM(argv[0]);
CHECK_COMPLEX(argv[1]);
CHECK_MATRIX_COMPLEX(argv[2]);
CHECK_VECTOR_COMPLEX(argv[3]);
type = FIX2INT(argv[0]);
Data_Get_Struct(argv[1], gsl_complex, a);
Data_Get_Struct(argv[2], gsl_matrix_complex, A);
Data_Get_Struct(argv[3], gsl_vector_complex, x);
istart = 4;
break;
default:
Data_Get_Struct(obj, gsl_matrix_complex, A);
if (argc < 3) rb_raise(rb_eArgError, "wrong number of arguments (%d for >= 3)",
argc);
CHECK_FIXNUM(argv[0]);
CHECK_COMPLEX(argv[1]);
CHECK_VECTOR_COMPLEX(argv[2]);
type = FIX2INT(argv[0]);
Data_Get_Struct(argv[1], gsl_complex, a);
Data_Get_Struct(argv[2], gsl_vector_complex, x);
istart = 3;
break;
}
switch (argc - istart) {
case 2:
CHECK_COMPLEX(argv[istart]);
CHECK_VECTOR_COMPLEX(argv[istart+1]);
Data_Get_Struct(argv[istart], gsl_complex, b);
Data_Get_Struct(argv[istart+1], gsl_vector_complex, y);
break;
case 0:
z = gsl_complex_rect(0.0, 0.0);
b = &z;
y = gsl_vector_complex_alloc(x->size);
flag = 1;
break;
default:
rb_raise(rb_eArgError, "wrong number of arguments");
break;
}
gsl_blas_zgemv(type, *a, A, x, *b, y);
if (flag == 1) return Data_Wrap_Struct(cgsl_vector_complex, 0, gsl_vector_complex_free, y);
else return argv[argc-1];
}
static VALUE rb_gsl_blas_zgemv2(int argc, VALUE *argv, VALUE obj)
{
gsl_matrix_complex *A = NULL;
gsl_vector_complex *x = NULL, *y, *ynew;
gsl_complex *a, *b, z;
int type, istart, flag = 0;
switch (TYPE(obj)) {
case T_MODULE:
case T_CLASS:
case T_OBJECT:
if (argc < 4) rb_raise(rb_eArgError, "wrong number of arguments (%d for >= 4)",
argc);
CHECK_FIXNUM(argv[0]);
CHECK_COMPLEX(argv[1]);
CHECK_MATRIX_COMPLEX(argv[2]);
CHECK_VECTOR_COMPLEX(argv[3]);
type = FIX2INT(argv[0]);
Data_Get_Struct(argv[1], gsl_complex, a);
Data_Get_Struct(argv[2], gsl_matrix_complex, A);
Data_Get_Struct(argv[3], gsl_vector_complex, x);
istart = 4;
break;
default:
Data_Get_Struct(obj, gsl_matrix_complex, A);
if (argc < 3) rb_raise(rb_eArgError, "wrong number of arguments (%d for >= 3)",
argc);
CHECK_FIXNUM(argv[0]);
CHECK_COMPLEX(argv[1]);
CHECK_VECTOR_COMPLEX(argv[2]);
type = FIX2INT(argv[0]);
Data_Get_Struct(argv[1], gsl_complex, a);
Data_Get_Struct(argv[2], gsl_vector_complex, x);
istart = 3;
break;
}
switch (argc - istart) {
case 2:
CHECK_COMPLEX(argv[istart]);
CHECK_VECTOR_COMPLEX(argv[istart+1]);
Data_Get_Struct(argv[istart], gsl_complex, b);
Data_Get_Struct(argv[istart+1], gsl_vector_complex, y);
break;
case 0:
z = gsl_complex_rect(0.0, 0.0);
b = &z;
y = gsl_vector_complex_alloc(x->size);
flag = 1;
break;
default:
rb_raise(rb_eArgError, "wrong number of arguments");
break;
}
ynew = gsl_vector_complex_alloc(y->size);
gsl_vector_complex_memcpy(ynew, y);
gsl_blas_zgemv(type, *a, A, x, *b, ynew);
if (flag == 1) gsl_vector_complex_free(y);
return Data_Wrap_Struct(cgsl_vector_complex, 0, gsl_vector_complex_free, ynew);
}
static VALUE rb_gsl_blas_dtrmv(int argc, VALUE *argv, VALUE obj)
{
gsl_matrix *A = NULL;
gsl_vector *x = NULL;
switch (TYPE(obj)) {
case T_MODULE:
case T_CLASS:
case T_OBJECT:
if (argc != 5) rb_raise(rb_eArgError, "wrong number of arguments (%d for 5)",
argc);
CHECK_MATRIX(argv[3]);
CHECK_VECTOR(argv[4]);
Data_Get_Struct(argv[3], gsl_matrix, A);
Data_Get_Struct(argv[4], gsl_vector, x);
break;
default:
if (argc != 4) rb_raise(rb_eArgError, "wrong number of arguments (%d for 4)",
argc);
CHECK_VECTOR(argv[3]);
Data_Get_Struct(obj, gsl_matrix, A);
Data_Get_Struct(argv[3], gsl_vector, x);
break;
}
CHECK_FIXNUM(argv[0]); CHECK_FIXNUM(argv[1]); CHECK_FIXNUM(argv[2]);
gsl_blas_dtrmv(FIX2INT(argv[0]), FIX2INT(argv[1]), FIX2INT(argv[2]),
A, x);
return argv[argc-1];
}
static VALUE rb_gsl_blas_dtrmv2(int argc, VALUE *argv, VALUE obj)
{
gsl_matrix *A = NULL;
gsl_vector *x = NULL, *xnew;
switch (TYPE(obj)) {
case T_MODULE:
case T_CLASS:
case T_OBJECT:
if (argc != 5) rb_raise(rb_eArgError, "wrong number of arguments (%d for 5)",
argc);
CHECK_MATRIX(argv[3]);
CHECK_VECTOR(argv[4]);
Data_Get_Struct(argv[3], gsl_matrix, A);
Data_Get_Struct(argv[4], gsl_vector, x);
break;
default:
if (argc != 4) rb_raise(rb_eArgError, "wrong number of arguments (%d for 4)",
argc);
CHECK_VECTOR(argv[3]);
Data_Get_Struct(obj, gsl_matrix, A);
Data_Get_Struct(argv[3], gsl_vector, x);
break;
}
CHECK_FIXNUM(argv[0]); CHECK_FIXNUM(argv[1]); CHECK_FIXNUM(argv[2]);
xnew = gsl_vector_alloc(x->size);
gsl_vector_memcpy(xnew, x);
gsl_blas_dtrmv(FIX2INT(argv[0]), FIX2INT(argv[1]), FIX2INT(argv[2]),
A, xnew);
return Data_Wrap_Struct(cgsl_vector, 0, gsl_vector_free, xnew);
}
static VALUE rb_gsl_blas_ztrmv(int argc, VALUE *argv, VALUE obj)
{
gsl_matrix_complex *A = NULL;
gsl_vector_complex *x = NULL;
switch (TYPE(obj)) {
case T_MODULE:
case T_CLASS:
case T_OBJECT:
if (argc != 5) rb_raise(rb_eArgError, "wrong number of arguments (%d for 5)",
argc);
CHECK_MATRIX_COMPLEX(argv[3]);
CHECK_VECTOR_COMPLEX(argv[4]);
Data_Get_Struct(argv[3], gsl_matrix_complex, A);
Data_Get_Struct(argv[4], gsl_vector_complex, x);
break;
default:
if (argc != 4) rb_raise(rb_eArgError, "wrong number of arguments (%d for 4)",
argc);
CHECK_VECTOR_COMPLEX(argv[3]);
Data_Get_Struct(obj, gsl_matrix_complex, A);
Data_Get_Struct(argv[3], gsl_vector_complex, x);
break;
}
CHECK_FIXNUM(argv[0]); CHECK_FIXNUM(argv[1]); CHECK_FIXNUM(argv[2]);
gsl_blas_ztrmv(FIX2INT(argv[0]), FIX2INT(argv[1]), FIX2INT(argv[2]),
A, x);
return argv[argc-1];
}
static VALUE rb_gsl_blas_ztrmv2(int argc, VALUE *argv, VALUE obj)
{
gsl_matrix_complex *A = NULL;
gsl_vector_complex *x = NULL, *xnew;
switch (TYPE(obj)) {
case T_MODULE:
case T_CLASS:
case T_OBJECT:
if (argc != 5) rb_raise(rb_eArgError, "wrong number of arguments (%d for 5)",
argc);
CHECK_MATRIX_COMPLEX(argv[3]);
CHECK_VECTOR_COMPLEX(argv[4]);
Data_Get_Struct(argv[3], gsl_matrix_complex, A);
Data_Get_Struct(argv[4], gsl_vector_complex, x);
break;
default:
if (argc != 4) rb_raise(rb_eArgError, "wrong number of arguments (%d for 4)",
argc);
CHECK_VECTOR_COMPLEX(argv[3]);
Data_Get_Struct(obj, gsl_matrix_complex, A);
Data_Get_Struct(argv[3], gsl_vector_complex, x);
break;
}
CHECK_FIXNUM(argv[0]); CHECK_FIXNUM(argv[1]); CHECK_FIXNUM(argv[2]);
xnew = gsl_vector_complex_alloc(x->size);
gsl_vector_complex_memcpy(xnew, x);
gsl_blas_ztrmv(FIX2INT(argv[0]), FIX2INT(argv[1]), FIX2INT(argv[2]),
A, xnew);
return Data_Wrap_Struct(cgsl_vector_complex, 0, gsl_vector_complex_free, xnew);
}
static VALUE rb_gsl_blas_dtrsv(int argc, VALUE *argv, VALUE obj)
{
gsl_matrix *A = NULL;
gsl_vector *x = NULL;
switch (TYPE(obj)) {
case T_MODULE:
case T_CLASS:
case T_OBJECT:
if (argc != 5) rb_raise(rb_eArgError, "wrong number of arguments (%d for 5)",
argc);
CHECK_MATRIX(argv[3]);
CHECK_VECTOR(argv[4]);
Data_Get_Struct(argv[3], gsl_matrix, A);
Data_Get_Struct(argv[4], gsl_vector, x);
break;
default:
if (argc != 4) rb_raise(rb_eArgError, "wrong number of arguments (%d for 4)",
argc);
CHECK_VECTOR(argv[3]);
Data_Get_Struct(obj, gsl_matrix, A);
Data_Get_Struct(argv[3], gsl_vector, x);
break;
}
CHECK_FIXNUM(argv[0]); CHECK_FIXNUM(argv[1]); CHECK_FIXNUM(argv[2]);
gsl_blas_dtrsv(FIX2INT(argv[0]), FIX2INT(argv[1]), FIX2INT(argv[2]),
A, x);
return argv[argc-1];
}
static VALUE rb_gsl_blas_dtrsv2(int argc, VALUE *argv, VALUE obj)
{
gsl_matrix *A = NULL;
gsl_vector *x = NULL, *xnew;
switch (TYPE(obj)) {
case T_MODULE:
case T_CLASS:
case T_OBJECT:
if (argc != 5) rb_raise(rb_eArgError, "wrong number of arguments (%d for 5)",
argc);
CHECK_MATRIX(argv[3]);
CHECK_VECTOR(argv[4]);
Data_Get_Struct(argv[3], gsl_matrix, A);
Data_Get_Struct(argv[4], gsl_vector, x);
break;
default:
if (argc != 4) rb_raise(rb_eArgError, "wrong number of arguments (%d for 4)",
argc);
CHECK_VECTOR(argv[3]);
Data_Get_Struct(obj, gsl_matrix, A);
Data_Get_Struct(argv[3], gsl_vector, x);
break;
}
CHECK_FIXNUM(argv[0]); CHECK_FIXNUM(argv[1]); CHECK_FIXNUM(argv[2]);
xnew = gsl_vector_alloc(x->size);
gsl_vector_memcpy(xnew, x);
gsl_blas_dtrsv(FIX2INT(argv[0]), FIX2INT(argv[1]), FIX2INT(argv[2]),
A, xnew);
return Data_Wrap_Struct(cgsl_vector, 0, gsl_vector_free, xnew);
}
static VALUE rb_gsl_blas_ztrsv(int argc, VALUE *argv, VALUE obj)
{
gsl_matrix_complex *A = NULL;
gsl_vector_complex *x = NULL;
switch (TYPE(obj)) {
case T_MODULE:
case T_CLASS:
case T_OBJECT:
if (argc != 5) rb_raise(rb_eArgError, "wrong number of arguments (%d for 5)",
argc);
CHECK_MATRIX_COMPLEX(argv[3]);
CHECK_VECTOR_COMPLEX(argv[4]);
Data_Get_Struct(argv[3], gsl_matrix_complex, A);
Data_Get_Struct(argv[4], gsl_vector_complex, x);
break;
default:
if (argc != 4) rb_raise(rb_eArgError, "wrong number of arguments (%d for 4)",
argc);
CHECK_VECTOR_COMPLEX(argv[3]);
Data_Get_Struct(obj, gsl_matrix_complex, A);
Data_Get_Struct(argv[3], gsl_vector_complex, x);
break;
}
CHECK_FIXNUM(argv[0]); CHECK_FIXNUM(argv[1]); CHECK_FIXNUM(argv[2]);
gsl_blas_ztrsv(FIX2INT(argv[0]), FIX2INT(argv[1]), FIX2INT(argv[2]),
A, x);
return argv[argc-1];
}
static VALUE rb_gsl_blas_ztrsv2(int argc, VALUE *argv, VALUE obj)
{
gsl_matrix_complex *A = NULL;
gsl_vector_complex *x = NULL, *xnew;
switch (TYPE(obj)) {
case T_MODULE:
case T_CLASS:
case T_OBJECT:
if (argc != 5) rb_raise(rb_eArgError, "wrong number of arguments (%d for 5)",
argc);
CHECK_MATRIX_COMPLEX(argv[3]);
CHECK_VECTOR_COMPLEX(argv[4]);
Data_Get_Struct(argv[3], gsl_matrix_complex, A);
Data_Get_Struct(argv[4], gsl_vector_complex, x);
break;
default:
if (argc != 4) rb_raise(rb_eArgError, "wrong number of arguments (%d for 4)",
argc);
CHECK_VECTOR_COMPLEX(argv[3]);
Data_Get_Struct(obj, gsl_matrix_complex, A);
Data_Get_Struct(argv[3], gsl_vector_complex, x);
break;
}
CHECK_FIXNUM(argv[0]); CHECK_FIXNUM(argv[1]); CHECK_FIXNUM(argv[2]);
xnew = gsl_vector_complex_alloc(x->size);
gsl_vector_complex_memcpy(xnew, x);
gsl_blas_ztrsv(FIX2INT(argv[0]), FIX2INT(argv[1]), FIX2INT(argv[2]),
A, xnew);
return Data_Wrap_Struct(cgsl_vector_complex, 0, gsl_vector_complex_free, xnew);
}
static VALUE rb_gsl_blas_dsymv(int argc, VALUE *argv, VALUE obj)
{
gsl_matrix *A = NULL;
gsl_vector *x = NULL, *y = NULL;
double a, b;
int type, istart, flag = 0;
switch (TYPE(obj)) {
case T_MODULE:
case T_CLASS:
case T_OBJECT:
if (argc < 4) rb_raise(rb_eArgError, "wrong number of arguments (%d for >= 4)",
argc);
CHECK_FIXNUM(argv[0]);
Need_Float(argv[1]);
CHECK_MATRIX(argv[2]);
CHECK_VECTOR(argv[3]);
type = FIX2INT(argv[0]);
a = NUM2DBL(argv[1]);
Data_Get_Struct(argv[2], gsl_matrix, A);
Data_Get_Struct(argv[3], gsl_vector, x);
istart = 4;
break;
default:
Data_Get_Struct(obj, gsl_matrix, A);
if (argc < 3) rb_raise(rb_eArgError, "wrong number of arguments (%d for >= 3)",
argc);
CHECK_FIXNUM(argv[0]);
Need_Float(argv[1]);
CHECK_VECTOR(argv[2]);
type = FIX2INT(argv[0]);
a = NUM2DBL(argv[1]);
Data_Get_Struct(argv[2], gsl_vector, x);
istart = 3;
break;
}
switch (argc - istart) {
case 2:
Need_Float(argv[istart]);
CHECK_VECTOR(argv[istart+1]);
b = NUM2DBL(argv[istart]);
Data_Get_Struct(argv[istart+1], gsl_vector, y);
break;
case 0:
b = 0.0;
y = gsl_vector_alloc(x->size);
flag = 1;
break;
default:
rb_raise(rb_eArgError, "wrong number of arguments");
break;
}
gsl_blas_dsymv(type, a, A, x, b, y);
if (flag == 1) return Data_Wrap_Struct(cgsl_vector, 0, gsl_vector_free, y);
else return argv[argc-1];
}
static VALUE rb_gsl_blas_dsymv2(int argc, VALUE *argv, VALUE obj)
{
gsl_matrix *A = NULL;
gsl_vector *x = NULL, *y, *ynew;
double a, b;
int type, istart, flag = 0;
switch (TYPE(obj)) {
case T_MODULE:
case T_CLASS:
case T_OBJECT:
if (argc < 4) rb_raise(rb_eArgError, "wrong number of arguments (%d for >= 4)",
argc);
CHECK_FIXNUM(argv[0]);
Need_Float(argv[1]);
CHECK_MATRIX(argv[2]);
CHECK_VECTOR(argv[3]);
type = FIX2INT(argv[0]);
a = NUM2DBL(argv[1]);
Data_Get_Struct(argv[2], gsl_matrix, A);
Data_Get_Struct(argv[3], gsl_vector, x);
istart = 4;
break;
default:
Data_Get_Struct(obj, gsl_matrix, A);
if (argc < 3) rb_raise(rb_eArgError, "wrong number of arguments (%d for >= 3)",
argc);
CHECK_FIXNUM(argv[0]);
Need_Float(argv[1]);
CHECK_VECTOR(argv[2]);
type = FIX2INT(argv[0]);
a = NUM2DBL(argv[1]);
Data_Get_Struct(argv[2], gsl_vector, x);
istart = 3;
break;
}
switch (argc - istart) {
case 2:
Need_Float(argv[istart]);
CHECK_VECTOR(argv[istart+1]);
b = NUM2DBL(argv[istart]);
Data_Get_Struct(argv[istart+1], gsl_vector, y);
break;
case 0:
b = 0.0;
y = gsl_vector_alloc(x->size);
flag = 1;
break;
default:
rb_raise(rb_eArgError, "wrong number of arguments");
break;
}
ynew = gsl_vector_alloc(y->size);
gsl_vector_memcpy(ynew, y);
gsl_blas_dsymv(type, a, A, x, b, ynew);
if (flag == 1) gsl_vector_free(y);
return Data_Wrap_Struct(cgsl_vector, 0, gsl_vector_free, ynew);
}
static VALUE rb_gsl_blas_zhemv(int argc, VALUE *argv, VALUE obj)
{
gsl_matrix_complex *A = NULL;
gsl_vector_complex *x = NULL, *y = NULL;
gsl_complex *a, *b, z;
int type, istart, flag = 0;
switch (TYPE(obj)) {
case T_MODULE:
case T_CLASS:
case T_OBJECT:
if (argc < 4) rb_raise(rb_eArgError, "wrong number of arguments (%d for >= 4)",
argc);
CHECK_FIXNUM(argv[0]);
CHECK_COMPLEX(argv[1]);
CHECK_MATRIX_COMPLEX(argv[2]);
CHECK_VECTOR_COMPLEX(argv[3]);
type = FIX2INT(argv[0]);
Data_Get_Struct(argv[1], gsl_complex, a);
Data_Get_Struct(argv[2], gsl_matrix_complex, A);
Data_Get_Struct(argv[3], gsl_vector_complex, x);
istart = 4;
break;
default:
Data_Get_Struct(obj, gsl_matrix_complex, A);
if (argc < 3) rb_raise(rb_eArgError, "wrong number of arguments (%d for >= 3)",
argc);
CHECK_FIXNUM(argv[0]);
CHECK_COMPLEX(argv[1]);
CHECK_VECTOR_COMPLEX(argv[2]);
type = FIX2INT(argv[0]);
Data_Get_Struct(argv[1], gsl_complex, a);
Data_Get_Struct(argv[2], gsl_vector_complex, x);
istart = 3;
break;
}
switch (argc - istart) {
case 2:
CHECK_COMPLEX(argv[istart]);
CHECK_VECTOR_COMPLEX(argv[istart+1]);
Data_Get_Struct(argv[istart], gsl_complex, b);
Data_Get_Struct(argv[istart+1], gsl_vector_complex, y);
break;
case 0:
z = gsl_complex_rect(0.0, 0.0);
b = &z;
y = gsl_vector_complex_alloc(x->size);
flag = 1;
break;
default:
rb_raise(rb_eArgError, "wrong number of arguments");
break;
}
gsl_blas_zhemv(type, *a, A, x, *b, y);
if (flag == 1) return Data_Wrap_Struct(cgsl_vector_complex, 0, gsl_vector_complex_free, y);
else return argv[argc-1];
}
static VALUE rb_gsl_blas_zhemv2(int argc, VALUE *argv, VALUE obj)
{
gsl_matrix_complex *A = NULL;
gsl_vector_complex *x = NULL, *y, *ynew;
gsl_complex *a, *b, z;
int type, istart, flag = 0;
switch (TYPE(obj)) {
case T_MODULE:
case T_CLASS:
case T_OBJECT:
if (argc < 4) rb_raise(rb_eArgError, "wrong number of arguments (%d for >= 4)",
argc);
CHECK_FIXNUM(argv[0]);
CHECK_COMPLEX(argv[1]);
CHECK_MATRIX_COMPLEX(argv[2]);
CHECK_VECTOR_COMPLEX(argv[3]);
type = FIX2INT(argv[0]);
Data_Get_Struct(argv[1], gsl_complex, a);
Data_Get_Struct(argv[2], gsl_matrix_complex, A);
Data_Get_Struct(argv[3], gsl_vector_complex, x);
istart = 4;
break;
default:
Data_Get_Struct(obj, gsl_matrix_complex, A);
if (argc < 3) rb_raise(rb_eArgError, "wrong number of arguments (%d for >= 3)",
argc);
CHECK_FIXNUM(argv[0]);
CHECK_COMPLEX(argv[1]);
CHECK_VECTOR_COMPLEX(argv[2]);
type = FIX2INT(argv[0]);
Data_Get_Struct(argv[1], gsl_complex, a);
Data_Get_Struct(argv[2], gsl_vector_complex, x);
istart = 3;
break;
}
switch (argc - istart) {
case 2:
CHECK_COMPLEX(argv[istart]);
CHECK_VECTOR_COMPLEX(argv[istart+1]);
Data_Get_Struct(argv[istart], gsl_complex, b);
Data_Get_Struct(argv[istart+1], gsl_vector_complex, y);
break;
case 0:
z = gsl_complex_rect(0.0, 0.0);
b = &z;
y = gsl_vector_complex_alloc(x->size);
flag = 1;
break;
default:
rb_raise(rb_eArgError, "wrong number of arguments");
break;
}
ynew = gsl_vector_complex_alloc(y->size);
gsl_vector_complex_memcpy(ynew, y);
gsl_blas_zhemv(type, *a, A, x, *b, ynew);
if (flag == 1) gsl_vector_complex_free(y);
return Data_Wrap_Struct(cgsl_vector_complex, 0, gsl_vector_complex_free, ynew);
}
static VALUE rb_gsl_blas_dger(VALUE obj, VALUE aa, VALUE xx, VALUE yy, VALUE AA)
{
gsl_matrix *A = NULL;
gsl_vector *x = NULL, *y = NULL;
double a;
Need_Float(aa);
CHECK_VECTOR(xx); CHECK_VECTOR(yy);
CHECK_MATRIX(AA);
a = NUM2DBL(aa);
Data_Get_Struct(xx, gsl_vector, x);
Data_Get_Struct(yy, gsl_vector, y);
Data_Get_Struct(AA, gsl_matrix, A);
gsl_blas_dger(a, x, y, A);
return AA;
}
static VALUE rb_gsl_blas_dger2(VALUE obj, VALUE aa, VALUE xx, VALUE yy, VALUE AA)
{
gsl_matrix *A = NULL, *Anew = NULL;
gsl_vector *x = NULL, *y = NULL;
double a;
Need_Float(aa);
CHECK_VECTOR(xx); CHECK_VECTOR(yy);
CHECK_MATRIX(AA);
a = NUM2DBL(aa);
Data_Get_Struct(xx, gsl_vector, x);
Data_Get_Struct(yy, gsl_vector, y);
Data_Get_Struct(AA, gsl_matrix, A);
Anew = gsl_matrix_alloc(A->size1, A->size2);
gsl_matrix_memcpy(Anew, A);
gsl_blas_dger(a, x, y, Anew);
return Data_Wrap_Struct(cgsl_matrix, 0, gsl_matrix_free, Anew);
}
static VALUE rb_gsl_blas_zgeru(VALUE obj, VALUE aa, VALUE xx, VALUE yy, VALUE AA)
{
gsl_matrix_complex *A = NULL;
gsl_vector_complex *x = NULL, *y = NULL;
gsl_complex *a;
CHECK_COMPLEX(aa);
CHECK_VECTOR_COMPLEX(xx);CHECK_VECTOR_COMPLEX(yy);
CHECK_MATRIX_COMPLEX(AA);
Data_Get_Struct(aa, gsl_complex, a);
Data_Get_Struct(xx, gsl_vector_complex, x);
Data_Get_Struct(yy, gsl_vector_complex, y);
Data_Get_Struct(AA, gsl_matrix_complex, A);
gsl_blas_zgeru(*a, x, y, A);
return AA;
}
static VALUE rb_gsl_blas_zgeru2(VALUE obj, VALUE aa, VALUE xx, VALUE yy, VALUE AA)
{
gsl_matrix_complex *A = NULL, *Anew = NULL;
gsl_vector_complex *x = NULL, *y = NULL;
gsl_complex *a;
CHECK_COMPLEX(aa);
CHECK_VECTOR_COMPLEX(xx);CHECK_VECTOR_COMPLEX(yy);
CHECK_MATRIX_COMPLEX(AA);
Data_Get_Struct(aa, gsl_complex, a);
Data_Get_Struct(xx, gsl_vector_complex, x);
Data_Get_Struct(yy, gsl_vector_complex, y);
Data_Get_Struct(AA, gsl_matrix_complex, A);
Anew = gsl_matrix_complex_alloc(A->size1, A->size2);
gsl_matrix_complex_memcpy(Anew, A);
gsl_blas_zgeru(*a, x, y, Anew);
return Data_Wrap_Struct(cgsl_matrix_complex, 0, gsl_matrix_complex_free, Anew);
}
static VALUE rb_gsl_blas_zgerc(VALUE obj, VALUE aa, VALUE xx, VALUE yy, VALUE AA)
{
gsl_matrix_complex *A = NULL;
gsl_vector_complex *x = NULL, *y = NULL;
gsl_complex *a;
CHECK_COMPLEX(aa);
CHECK_VECTOR_COMPLEX(xx);CHECK_VECTOR_COMPLEX(yy);
CHECK_MATRIX_COMPLEX(AA);
Data_Get_Struct(aa, gsl_complex, a);
Data_Get_Struct(xx, gsl_vector_complex, x);
Data_Get_Struct(yy, gsl_vector_complex, y);
Data_Get_Struct(AA, gsl_matrix_complex, A);
gsl_blas_zgerc(*a, x, y, A);
return AA;
}
static VALUE rb_gsl_blas_zgerc2(VALUE obj, VALUE aa, VALUE xx, VALUE yy, VALUE AA)
{
gsl_matrix_complex *A = NULL, *Anew = NULL;
gsl_vector_complex *x = NULL, *y = NULL;
gsl_complex *a;
CHECK_COMPLEX(aa);
CHECK_VECTOR_COMPLEX(xx);CHECK_VECTOR_COMPLEX(yy);
CHECK_MATRIX_COMPLEX(AA);
Data_Get_Struct(aa, gsl_complex, a);
Data_Get_Struct(xx, gsl_vector_complex, x);
Data_Get_Struct(yy, gsl_vector_complex, y);
Data_Get_Struct(AA, gsl_matrix_complex, A);
Anew = gsl_matrix_complex_alloc(A->size1, A->size2);
gsl_matrix_complex_memcpy(Anew, A);
gsl_blas_zgerc(*a, x, y, Anew);
return Data_Wrap_Struct(cgsl_matrix_complex, 0, gsl_matrix_complex_free, Anew);
}
static VALUE rb_gsl_blas_dsyr(VALUE obj, VALUE tt, VALUE aa, VALUE xx, VALUE AA)
{
gsl_matrix *A = NULL;
gsl_vector *x = NULL;
double a;
CHECK_FIXNUM(tt);
Need_Float(aa);
CHECK_VECTOR(xx); CHECK_MATRIX(AA);
a = NUM2DBL(aa);
Data_Get_Struct(xx, gsl_vector, x);
Data_Get_Struct(AA, gsl_matrix, A);
gsl_blas_dsyr(FIX2INT(tt), a, x, A);
return AA;
}
static VALUE rb_gsl_blas_dsyr_a(VALUE obj, VALUE tt, VALUE aa, VALUE xx, VALUE AA)
{
gsl_matrix *A = NULL, *Anew = NULL;
gsl_vector *x = NULL;
double a;
CHECK_FIXNUM(tt);
Need_Float(aa);
CHECK_VECTOR(xx); CHECK_MATRIX(AA);
a = NUM2DBL(aa);
Data_Get_Struct(xx, gsl_vector, x);
Data_Get_Struct(AA, gsl_matrix, A);
Anew = gsl_matrix_alloc(A->size1, A->size2);
gsl_matrix_memcpy(Anew, A);
gsl_blas_dsyr(FIX2INT(tt), a, x, Anew);
return Data_Wrap_Struct(cgsl_matrix, 0, gsl_matrix_free, Anew);
}
static VALUE rb_gsl_blas_zher(VALUE obj, VALUE tt, VALUE aa, VALUE xx, VALUE AA)
{
gsl_matrix_complex *A = NULL;
gsl_vector_complex *x = NULL;
double a;
CHECK_FIXNUM(tt);
Need_Float(aa);
CHECK_VECTOR_COMPLEX(xx); CHECK_MATRIX_COMPLEX(AA);
a = NUM2DBL(aa);
Data_Get_Struct(xx, gsl_vector_complex, x);
Data_Get_Struct(AA, gsl_matrix_complex, A);
gsl_blas_zher(FIX2INT(tt), a, x, A);
return AA;
}
static VALUE rb_gsl_blas_zher_a(VALUE obj, VALUE tt, VALUE aa, VALUE xx, VALUE AA)
{
gsl_matrix_complex *A = NULL, *Anew = NULL;
gsl_vector_complex *x = NULL;
double a;
CHECK_FIXNUM(tt);
Need_Float(aa);
CHECK_VECTOR_COMPLEX(xx); CHECK_MATRIX_COMPLEX(AA);
a = NUM2DBL(aa);
Data_Get_Struct(xx, gsl_vector_complex, x);
Data_Get_Struct(AA, gsl_matrix_complex, A);
Anew = gsl_matrix_complex_alloc(A->size1, A->size2);
gsl_matrix_complex_memcpy(Anew, A);
gsl_blas_zher(FIX2INT(tt), a, x, Anew);
return Data_Wrap_Struct(cgsl_matrix_complex, 0, gsl_matrix_complex_free, Anew);
}
static VALUE rb_gsl_blas_dsyr2(VALUE obj, VALUE tt, VALUE aa, VALUE xx,
VALUE yy, VALUE AA)
{
gsl_matrix *A = NULL;
gsl_vector *x = NULL, *y = NULL;
double a;
CHECK_FIXNUM(tt);
Need_Float(aa);
CHECK_VECTOR(xx); CHECK_VECTOR(yy); CHECK_MATRIX(AA);
a = NUM2DBL(aa);
Data_Get_Struct(xx, gsl_vector, x);
Data_Get_Struct(yy, gsl_vector, y);
Data_Get_Struct(AA, gsl_matrix, A);
gsl_blas_dsyr2(FIX2INT(tt), a, x, y, A);
return AA;
}
static VALUE rb_gsl_blas_dsyr2_a(VALUE obj, VALUE tt, VALUE aa, VALUE xx,
VALUE yy, VALUE AA)
{
gsl_matrix *A = NULL, *Anew = NULL;
gsl_vector *x = NULL, *y = NULL;
double a;
CHECK_FIXNUM(tt);
Need_Float(aa);
CHECK_VECTOR(xx); CHECK_VECTOR(yy); CHECK_MATRIX(AA);
a = NUM2DBL(aa);
Data_Get_Struct(xx, gsl_vector, x);
Data_Get_Struct(yy, gsl_vector, y);
Data_Get_Struct(AA, gsl_matrix, A);
Anew = gsl_matrix_alloc(A->size1, A->size2);
gsl_matrix_memcpy(Anew, A);
gsl_blas_dsyr2(FIX2INT(tt), a, x, y, Anew);
return Data_Wrap_Struct(cgsl_matrix, 0, gsl_matrix_free, Anew);
}
static VALUE rb_gsl_blas_zher2(VALUE obj, VALUE tt, VALUE aa, VALUE xx,
VALUE yy, VALUE AA)
{
gsl_matrix_complex *A = NULL;
gsl_vector_complex *x = NULL, *y = NULL;
gsl_complex *a;
CHECK_FIXNUM(tt);
CHECK_COMPLEX(aa);
CHECK_VECTOR_COMPLEX(xx); CHECK_VECTOR_COMPLEX(yy); CHECK_MATRIX_COMPLEX(AA);
Data_Get_Struct(aa, gsl_complex, a);
Data_Get_Struct(xx, gsl_vector_complex, x);
Data_Get_Struct(yy, gsl_vector_complex, y);
Data_Get_Struct(AA, gsl_matrix_complex, A);
gsl_blas_zher2(FIX2INT(tt), *a, x, y, A);
return AA;
}
static VALUE rb_gsl_blas_zher2_a(VALUE obj, VALUE tt, VALUE aa, VALUE xx,
VALUE yy, VALUE AA)
{
gsl_matrix_complex *A = NULL, *Anew = NULL;
gsl_vector_complex *x = NULL, *y = NULL;
gsl_complex *a;
CHECK_FIXNUM(tt);
CHECK_COMPLEX(aa);
CHECK_VECTOR_COMPLEX(xx); CHECK_VECTOR_COMPLEX(yy); CHECK_MATRIX_COMPLEX(AA);
Data_Get_Struct(aa, gsl_complex, a);
Data_Get_Struct(xx, gsl_vector_complex, x);
Data_Get_Struct(yy, gsl_vector_complex, y);
Data_Get_Struct(AA, gsl_matrix_complex, A);
Anew = gsl_matrix_complex_alloc(A->size1, A->size2);
gsl_matrix_complex_memcpy(Anew, A);
gsl_blas_zher2(FIX2INT(tt), *a, x, y, Anew);
return Data_Wrap_Struct(cgsl_matrix_complex, 0, gsl_matrix_complex_free, Anew);
}
void Init_gsl_blas2(VALUE module)
{
define_const(module);
rb_define_module_function(module, "dgemv!", rb_gsl_blas_dgemv, -1);
rb_define_method(cgsl_matrix, "blas_dgemv!", rb_gsl_blas_dgemv, -1);
rb_define_alias(cgsl_matrix, "dgemv!", "blas_dgemv!");
rb_define_alias(cgsl_matrix, "gemv!", "blas_dgemv!");
rb_define_module_function(module, "dgemv", rb_gsl_blas_dgemv2, -1);
rb_define_method(cgsl_matrix, "blas_dgemv", rb_gsl_blas_dgemv2, -1);
rb_define_alias(cgsl_matrix, "dgemv", "blas_dgemv");
rb_define_alias(cgsl_matrix, "gemv", "blas_dgemv");
rb_define_module_function(module, "zgemv!", rb_gsl_blas_zgemv, -1);
rb_define_method(cgsl_matrix_complex, "blas_zgemv!", rb_gsl_blas_zgemv, -1);
rb_define_alias(cgsl_matrix_complex, "zgemv!", "blas_zgemv!");
rb_define_alias(cgsl_matrix_complex, "gemv!", "blas_zgemv!");
rb_define_module_function(module, "zgemv", rb_gsl_blas_zgemv2, -1);
rb_define_method(cgsl_matrix_complex, "blas_zgemv", rb_gsl_blas_zgemv2, -1);
rb_define_alias(cgsl_matrix_complex, "zgemv", "blas_zgemv");
rb_define_alias(cgsl_matrix_complex, "gemv", "blas_zgemv");
rb_define_module_function(module, "dtrmv!", rb_gsl_blas_dtrmv, -1);
rb_define_method(cgsl_matrix, "blas_dtrmv!", rb_gsl_blas_dtrmv, -1);
rb_define_alias(cgsl_matrix, "dtrmv!", "blas_dtrmv!");
rb_define_alias(cgsl_matrix, "trmv!", "blas_dtrmv!");
rb_define_module_function(module, "dtrmv", rb_gsl_blas_dtrmv2, -1);
rb_define_method(cgsl_matrix, "blas_dtrmv", rb_gsl_blas_dtrmv2, -1);
rb_define_alias(cgsl_matrix, "dtrmv", "blas_dtrmv");
rb_define_alias(cgsl_matrix, "trmv", "blas_dtrmv");
rb_define_module_function(module, "ztrmv!", rb_gsl_blas_ztrmv, -1);
rb_define_method(cgsl_matrix_complex, "blas_ztrmv!", rb_gsl_blas_ztrmv, -1);
rb_define_alias(cgsl_matrix_complex, "ztrmv!", "blas_ztrmv!");
rb_define_module_function(module, "ztrmv", rb_gsl_blas_ztrmv2, -1);
rb_define_method(cgsl_matrix_complex, "blas_ztrmv", rb_gsl_blas_ztrmv2, -1);
rb_define_alias(cgsl_matrix_complex, "ztrmv", "blas_ztrmv");
rb_define_alias(cgsl_matrix_complex, "trmv", "blas_ztrmv");
rb_define_module_function(module, "dtrsv!", rb_gsl_blas_dtrsv, -1);
rb_define_method(cgsl_matrix, "blas_dtrsv!", rb_gsl_blas_dtrsv, -1);
rb_define_alias(cgsl_matrix, "dtrsv!", "blas_dtrsv!");
rb_define_alias(cgsl_matrix, "trsv!", "blas_dtrsv!");
rb_define_module_function(module, "dtrsv", rb_gsl_blas_dtrsv2, -1);
rb_define_method(cgsl_matrix, "blas_dtrsv", rb_gsl_blas_dtrsv2, -1);
rb_define_alias(cgsl_matrix, "dtrsv", "blas_dtrsv");
rb_define_alias(cgsl_matrix, "trsv", "blas_dtrsv");
rb_define_module_function(module, "ztrsv!", rb_gsl_blas_ztrsv, -1);
rb_define_method(cgsl_matrix_complex, "blas_ztrsv!", rb_gsl_blas_ztrsv, -1);
rb_define_alias(cgsl_matrix_complex, "ztrsv!", "blas_ztrsv!");
rb_define_alias(cgsl_matrix_complex, "trsv!", "blas_ztrsv!");
rb_define_module_function(module, "ztrsv", rb_gsl_blas_ztrsv2, -1);
rb_define_method(cgsl_matrix_complex, "blas_ztrsv", rb_gsl_blas_ztrsv2, -1);
rb_define_alias(cgsl_matrix_complex, "ztrsv", "blas_ztrsv");
rb_define_alias(cgsl_matrix_complex, "trsv", "blas_ztrsv");
rb_define_module_function(module, "dsymv!", rb_gsl_blas_dsymv, -1);
rb_define_method(cgsl_matrix, "blas_dsymv!", rb_gsl_blas_dsymv, -1);
rb_define_alias(cgsl_matrix, "dsymv!", "blas_dsymv!");
rb_define_alias(cgsl_matrix, "symv!", "blas_dsymv!");
rb_define_module_function(module, "dsymv", rb_gsl_blas_dsymv2, -1);
rb_define_method(cgsl_matrix, "blas_dsymv", rb_gsl_blas_dsymv2, -1);
rb_define_alias(cgsl_matrix, "dsymv", "blas_dsymv");
rb_define_alias(cgsl_matrix, "symv", "blas_dsymv");
rb_define_module_function(module, "zhemv!", rb_gsl_blas_zhemv, -1);
rb_define_method(cgsl_matrix_complex, "blas_zhemv!", rb_gsl_blas_zhemv, -1);
rb_define_alias(cgsl_matrix_complex, "zhemv!", "blas_zhemv!");
rb_define_alias(cgsl_matrix_complex, "symv!", "blas_zhemv!");
rb_define_module_function(module, "zhemv", rb_gsl_blas_zhemv2, -1);
rb_define_method(cgsl_matrix_complex, "blas_zhemv", rb_gsl_blas_zhemv2, -1);
rb_define_alias(cgsl_matrix_complex, "zhemv", "blas_zhemv");
rb_define_alias(cgsl_matrix_complex, "symv", "blas_zhemv");
rb_define_module_function(module, "dger!", rb_gsl_blas_dger, 4);
rb_define_module_function(module, "dger", rb_gsl_blas_dger2, 4);
rb_define_module_function(module, "zgeru!", rb_gsl_blas_zgeru, 4);
rb_define_module_function(module, "zgeru", rb_gsl_blas_zgeru2, 4);
rb_define_module_function(module, "zgerc!", rb_gsl_blas_zgerc, 4);
rb_define_module_function(module, "zgerc", rb_gsl_blas_zgerc2, 4);
rb_define_module_function(module, "dsyr!", rb_gsl_blas_dsyr, 4);
rb_define_module_function(module, "dsyr", rb_gsl_blas_dsyr_a, 4);
rb_define_module_function(module, "zher!", rb_gsl_blas_zher, 4);
rb_define_module_function(module, "zher", rb_gsl_blas_zher_a, 4);
rb_define_module_function(module, "dsyr2!", rb_gsl_blas_dsyr2, 4);
rb_define_module_function(module, "dsyr2", rb_gsl_blas_dsyr2_a, 4);
rb_define_module_function(module, "zher2!", rb_gsl_blas_zher2, 4);
rb_define_module_function(module, "zher2", rb_gsl_blas_zher2_a, 4);
}
syntax highlighted by Code2HTML, v. 0.9.1