C++ CHECK_MATRIX_COMPLEX函数代码示例

static VALUE rb_Dirac_matrix_is_equal(int argc, VALUE *argv, VALUE obj)
{
  gsl_complex ztmp, *z;
  gsl_matrix_complex *m1, *m2;
  VALUE vz;
  switch (TYPE(obj)) {
  case T_MODULE:
  case T_CLASS:
  case T_OBJECT:
    CHECK_MATRIX_COMPLEX(argv[0]);
    CHECK_MATRIX_COMPLEX(argv[1]);
    Data_Get_Struct(argv[0], gsl_matrix_complex, m1);
    Data_Get_Struct(argv[1], gsl_matrix_complex, m2);
    if (matrix_is_equal(m1, m2, &ztmp)) {
      vz = Data_Make_Struct(cgsl_complex, gsl_complex, 0, free, z);
      *z = ztmp;
      return vz;
    } else {
      return Qfalse;
    }
    break;
  default:
    CHECK_MATRIX_COMPLEX(argv[0]);
    Data_Get_Struct(obj, gsl_matrix_complex, m1);
    Data_Get_Struct(argv[0], gsl_matrix_complex, m2);
    if (matrix_is_equal(m1, m2, &ztmp)) {
      vz = Data_Make_Struct(cgsl_complex, gsl_complex, 0, free, z);
      *z = ztmp;
      return vz;
    } else {
      return Qfalse;
    }
    break;
  }
}

static VALUE rb_gsl_blas_ztrsm2(VALUE obj, VALUE s, VALUE u, VALUE ta,
             VALUE d, VALUE a, VALUE aa, VALUE bb)
{
  gsl_matrix_complex *A = NULL, *B = NULL, *Bnew = NULL;
  gsl_complex *pa = NULL;
  CBLAS_SIDE_t Side;
  CBLAS_UPLO_t Uplo;
  CBLAS_TRANSPOSE_t TransA;
  CBLAS_DIAG_t Diag;
  CHECK_FIXNUM(s);  CHECK_FIXNUM(u);  CHECK_FIXNUM(ta);  CHECK_FIXNUM(d);
  CHECK_COMPLEX(a);
  CHECK_MATRIX_COMPLEX(aa);
  CHECK_MATRIX_COMPLEX(bb);
  Side = FIX2INT(s);
  Uplo = FIX2INT(u);
  TransA = FIX2INT(ta);
  Diag = FIX2INT(d);
  Data_Get_Struct(a, gsl_complex, pa);
  Data_Get_Struct(aa, gsl_matrix_complex, A);
  Data_Get_Struct(bb, gsl_matrix_complex, B);
  Bnew = gsl_matrix_complex_alloc(B->size1, B->size2);
  gsl_matrix_complex_memcpy(Bnew, B);
  gsl_blas_ztrsm(Side, Uplo, TransA, Diag, *pa, A, Bnew);
  return Data_Wrap_Struct(cgsl_matrix_complex, 0, gsl_matrix_complex_free, Bnew);
}

static VALUE rb_gsl_linalg_complex_LU_refine(VALUE obj, VALUE vm,
					     VALUE lu, VALUE pp, VALUE bb,
					     VALUE xx)
{
  gsl_matrix_complex *m = NULL, *mlu = NULL;
  gsl_permutation *p = NULL;
  gsl_vector_complex *b = NULL, *x = NULL, *r = NULL;
  int flagb = 0;
  VALUE vr;

  if (CLASS_OF(obj) != cgsl_matrix_complex_LU)
    rb_raise(rb_eRuntimeError, "Decompose first!");
  CHECK_MATRIX_COMPLEX(vm);
  CHECK_MATRIX_COMPLEX(lu);
  CHECK_PERMUTATION(pp);
  CHECK_VECTOR_COMPLEX(xx);
  Data_Get_Struct(vm, gsl_matrix_complex, m);
  Data_Get_Struct(lu, gsl_matrix_complex, mlu);
  Data_Get_Struct(pp, gsl_permutation, p);
  CHECK_VECTOR_COMPLEX(bb);
  Data_Get_Struct(bb, gsl_vector_complex, b);
  Data_Get_Struct(xx, gsl_vector_complex, x);
  r = gsl_vector_complex_alloc(m->size1);
  gsl_linalg_complex_LU_refine(m, mlu, p, b, x, r);
  vr = Data_Wrap_Struct(cgsl_vector_complex, 0, gsl_vector_complex_free, r);
  if (flagb == 1) gsl_vector_complex_free(b);
  return rb_ary_new3(2, xx, vr);
}

static VALUE rb_gsl_linalg_cholesky_svx(int argc, VALUE *argv, VALUE obj)
{
  gsl_matrix_complex *A = NULL, *Atmp = NULL;
  gsl_vector_complex *b = NULL;
  int flaga = 0;
  VALUE vA, vb;
  switch(TYPE(obj)) {
  case T_MODULE:  case T_CLASS:  case T_OBJECT:
    if (argc != 2) rb_raise(rb_eArgError, "wrong number of argument (%d for 2)",
			    argc);
    vA = argv[0];
    vb = argv[1];
    break;
  default:
    if (argc != 1) rb_raise(rb_eArgError, "wrong number of argument (%d for 1)",
			    argc);
    vA = obj;
    vb = argv[0];
    break;
  }
  CHECK_MATRIX_COMPLEX(vA);
  Data_Get_Struct(vA, gsl_matrix_complex, Atmp);
  CHECK_VECTOR_COMPLEX(vb);
  Data_Get_Struct(vb, gsl_vector_complex, b);
  if (CLASS_OF(vA) == cgsl_matrix_complex_C) {
    A = Atmp;
  } else {
    A = make_matrix_complex_clone(Atmp);
    flaga = 1;
    gsl_linalg_complex_cholesky_decomp(A);
  }
  gsl_linalg_complex_cholesky_svx(A, b);
  if (flaga == 1) gsl_matrix_complex_free(A);
  return vb;
}

static VALUE rb_dirac_anticommute(VALUE obj, VALUE mm1, VALUE mm2)
{
  gsl_matrix_complex *m1, *m2;
  gsl_matrix_complex *mnew1, *mnew2;
  CHECK_MATRIX_COMPLEX(mm1);
  CHECK_MATRIX_COMPLEX(mm2);
  Data_Get_Struct(mm1, gsl_matrix_complex, m1);
  Data_Get_Struct(mm2, gsl_matrix_complex, m2);
  mnew1 = gsl_matrix_complex_alloc(m1->size1, m1->size2);
  mnew2 = gsl_matrix_complex_alloc(m1->size1, m1->size2);
  gsl_matrix_complex_mul(mnew1, m1, m2);
  gsl_matrix_complex_mul(mnew2, m2, m1);
  gsl_matrix_complex_add(mnew1, mnew2);
  gsl_matrix_complex_free(mnew2);
  return Data_Wrap_Struct(cgsl_matrix_complex, 0, gsl_matrix_complex_free,
                          mnew1);
}

VALUE rb_gsl_linalg_complex_LU_decomp(int argc, VALUE *argv, VALUE obj)
{
  gsl_matrix_complex *m = NULL;
  gsl_permutation *p = NULL;
  int signum, itmp;
  size_t size;
  VALUE obj2;

  switch (TYPE(obj)) {
  case T_MODULE:
  case T_CLASS:
  case T_OBJECT:
    if (argc != 1) rb_raise(rb_eArgError, "wrong number of arguments (%d for 1)",
			    argc);
    CHECK_MATRIX_COMPLEX(argv[0]);
    Data_Get_Struct(argv[0], gsl_matrix_complex, m);
    itmp = 1;
    break;
  default:
    CHECK_MATRIX_COMPLEX(obj);
    Data_Get_Struct(obj, gsl_matrix_complex, m);
    itmp = 0;
  }
  size = m->size1;
  switch (argc-itmp) {
  case 0:
    p = gsl_permutation_alloc(size);
    gsl_linalg_complex_LU_decomp(m, p, &signum);
    if (itmp == 1) rb_obj_reveal(argv[0], cgsl_matrix_complex_LU);
    else rb_obj_reveal(obj, cgsl_matrix_complex_LU);
    obj2 = Data_Wrap_Struct(cgsl_permutation, 0, gsl_permutation_free, p);
    return rb_ary_new3(2, obj2, INT2FIX(signum));
    break;
  case 1:  /* when a permutation object is given */
    CHECK_PERMUTATION(argv[itmp]);
    Data_Get_Struct(argv[itmp], gsl_permutation, p);
    gsl_linalg_complex_LU_decomp(m, p, &signum);
    if (itmp == 1) rb_obj_reveal(argv[0], cgsl_matrix_complex_LU);
    else rb_obj_reveal(obj, cgsl_matrix_complex_LU);
    return INT2FIX(signum);
    break;
  default:
    rb_raise(rb_eArgError, "Usage: LU_decomp!() or LU_decomp!(permutation)");
  }
}

static VALUE rb_gsl_blas_zsyrk(VALUE obj, VALUE u, VALUE t, VALUE a, VALUE aa,
             VALUE b, VALUE cc)
{
  gsl_matrix_complex *A = NULL, *C = NULL;
  gsl_complex *pa = NULL, *pb = NULL;
  CBLAS_UPLO_t Uplo;
  CBLAS_TRANSPOSE_t Trans;
  CHECK_FIXNUM(u);  CHECK_FIXNUM(t);
  CHECK_COMPLEX(a);  CHECK_COMPLEX(b);
  CHECK_MATRIX_COMPLEX(aa);  CHECK_MATRIX_COMPLEX(cc);
  Uplo = FIX2INT(u);
  Trans = FIX2INT(t);
  Data_Get_Struct(a, gsl_complex, pa);
  Data_Get_Struct(b, gsl_complex, pb);
  Data_Get_Struct(aa, gsl_matrix_complex, A);
  Data_Get_Struct(cc, gsl_matrix_complex, C);
  gsl_blas_zsyrk(Uplo, Trans, *pa, A, *pb, C);
  return cc;
}

static VALUE rb_gsl_linalg_cholesky_decomp(int argc, VALUE *argv, VALUE obj)
{
  gsl_matrix_complex *A = NULL, *Atmp = NULL;
  switch(TYPE(obj)) {
  case T_MODULE:  case T_CLASS:  case T_OBJECT:
    if (argc != 1) rb_raise(rb_eArgError, "wrong number of argument (%d for 1)",
			    argc);
    CHECK_MATRIX_COMPLEX(argv[0]);
    Data_Get_Struct(argv[0], gsl_matrix_complex, Atmp);
    break;
  default:
    CHECK_MATRIX_COMPLEX(obj);
    Data_Get_Struct(obj, gsl_matrix_complex, Atmp);
    break;
  }
  A = make_matrix_complex_clone(Atmp);
  gsl_linalg_complex_cholesky_decomp(A);
  return Data_Wrap_Struct(cgsl_matrix_complex_C, 0, gsl_matrix_complex_free, A);
}

static VALUE rb_gsl_blas_zherk(VALUE obj, VALUE u, VALUE t, VALUE a, VALUE aa,
             VALUE b, VALUE cc)
{
  gsl_matrix_complex *A = NULL, *C = NULL;
  double alpha, beta;
  CBLAS_UPLO_t Uplo;
  CBLAS_TRANSPOSE_t Trans;
  CHECK_FIXNUM(u);  CHECK_FIXNUM(t);
  Need_Float(a);  Need_Float(b);
  CHECK_MATRIX_COMPLEX(aa);
  CHECK_MATRIX_COMPLEX(cc);
  Uplo = FIX2INT(u);
  Trans = FIX2INT(t);
  alpha = NUM2DBL(a);
  beta = NUM2DBL(b);
  Data_Get_Struct(aa, gsl_matrix_complex, A);
  Data_Get_Struct(cc, gsl_matrix_complex, C);
  gsl_blas_zherk(Uplo, Trans, alpha, A, beta, C);
  return cc;
}

static VALUE rb_gsl_blas_zherk2(VALUE obj, VALUE u, VALUE t, VALUE a, VALUE aa,
             VALUE b, VALUE cc)
{
  gsl_matrix_complex *A = NULL, *C = NULL, *Cnew = NULL;
  double alpha, beta;
  CBLAS_UPLO_t Uplo;
  CBLAS_TRANSPOSE_t Trans;
  CHECK_FIXNUM(u);  CHECK_FIXNUM(t);
  Need_Float(a);  Need_Float(b);
  CHECK_MATRIX_COMPLEX(aa);  CHECK_MATRIX_COMPLEX(cc);
  Uplo = FIX2INT(u);
  Trans = FIX2INT(t);
  alpha = NUM2DBL(a);
  beta = NUM2DBL(b);
  Data_Get_Struct(aa, gsl_matrix_complex, A);
  Data_Get_Struct(cc, gsl_matrix_complex, C);
  Cnew = gsl_matrix_complex_alloc(C->size1, C->size2);
  gsl_matrix_complex_memcpy(Cnew, C);
  gsl_blas_zherk(Uplo, Trans, alpha, A, beta, Cnew);
  return Data_Wrap_Struct(cgsl_matrix_complex, 0, gsl_matrix_complex_free, Cnew);
}

static VALUE rb_gsl_blas_zsyrk2(VALUE obj, VALUE u, VALUE t, VALUE a, VALUE aa,
             VALUE b, VALUE cc)
{
  gsl_matrix_complex *A = NULL, *C = NULL, *Cnew = NULL;
  gsl_complex *pa = NULL, *pb = NULL;
  CBLAS_UPLO_t Uplo;
  CBLAS_TRANSPOSE_t Trans;
  CHECK_FIXNUM(u);  CHECK_FIXNUM(t);
  CHECK_COMPLEX(a);  CHECK_COMPLEX(b);
  CHECK_MATRIX_COMPLEX(aa);  CHECK_MATRIX_COMPLEX(cc);
  Uplo = FIX2INT(u);
  Trans = FIX2INT(t);
  Data_Get_Struct(a, gsl_complex, pa);
  Data_Get_Struct(b, gsl_complex, pb);
  Data_Get_Struct(aa, gsl_matrix_complex, A);
  Data_Get_Struct(cc, gsl_matrix_complex, C);
  Cnew = gsl_matrix_complex_alloc(C->size1, C->size2);
  gsl_matrix_complex_memcpy(Cnew, C);
  gsl_blas_zsyrk(Uplo, Trans, *pa, A, *pb, Cnew);
  return Data_Wrap_Struct(cgsl_matrix_complex, 0, gsl_matrix_complex_free, Cnew);
}

static VALUE rb_gsl_blas_ztrsm(VALUE obj, VALUE s, VALUE u, VALUE ta,
             VALUE d, VALUE a, VALUE aa, VALUE bb)
{
  gsl_matrix_complex *A = NULL, *B = NULL;
  gsl_complex *pa = NULL;
  CBLAS_SIDE_t Side;
  CBLAS_UPLO_t Uplo;
  CBLAS_TRANSPOSE_t TransA;
  CBLAS_DIAG_t Diag;
  CHECK_FIXNUM(s);  CHECK_FIXNUM(u);  CHECK_FIXNUM(ta);  CHECK_FIXNUM(d);
  CHECK_COMPLEX(a);
  CHECK_MATRIX_COMPLEX(aa);
  CHECK_MATRIX_COMPLEX(bb);
  Side = FIX2INT(s);
  Uplo = FIX2INT(u);
  TransA = FIX2INT(ta);
  Diag = FIX2INT(d);
  Data_Get_Struct(a, gsl_complex, pa);
  Data_Get_Struct(aa, gsl_matrix_complex, A);
  Data_Get_Struct(bb, gsl_matrix_complex, B);
  gsl_blas_ztrsm(Side, Uplo, TransA, Diag, *pa, A, B);
  return bb;
}

static VALUE rb_gsl_linalg_complex_householder_mh(VALUE obj, VALUE t, VALUE vv, VALUE aa)
{
  gsl_vector_complex *v = NULL;
  gsl_complex *tau;
  gsl_matrix_complex *A = NULL;
  CHECK_COMPLEX(t);  
  CHECK_VECTOR_COMPLEX(vv);
  CHECK_MATRIX_COMPLEX(aa);
  Data_Get_Struct(t, gsl_complex, tau);
  Data_Get_Struct(vv, gsl_vector_complex, v);
  Data_Get_Struct(aa, gsl_matrix_complex, A);
  gsl_linalg_complex_householder_hm(*tau, v, A);
  return aa;
}

static VALUE rb_gsl_linalg_complex_LU_sgndet(int argc, VALUE *argv, VALUE obj)
{
  gsl_matrix_complex *m = NULL, *mtmp = NULL;
  gsl_permutation *p = NULL;
  gsl_complex *z = NULL;
  VALUE vz;
  int flagm = 0, signum, itmp;
  switch (TYPE(obj)) {
  case T_MODULE:
  case T_CLASS:
  case T_OBJECT:
    CHECK_MATRIX_COMPLEX(argv[0]);
    Data_Get_Struct(argv[0], gsl_matrix_complex, m);
    if (CLASS_OF(argv[0]) != cgsl_matrix_complex_LU) {
      mtmp = gsl_matrix_complex_alloc(m->size1, m->size2);
      gsl_matrix_complex_memcpy(mtmp, m);
      flagm = 1;
    } else {
      mtmp = m;
    }
    itmp = 1;
    break;
  default:
    Data_Get_Struct(obj, gsl_matrix_complex, m);
    if (CLASS_OF(obj) != cgsl_matrix_complex_LU) {
      mtmp = gsl_matrix_complex_alloc(m->size1, m->size2);
      gsl_matrix_complex_memcpy(mtmp, m);
      flagm = 1;
    } else {
      mtmp = m;
    }
    itmp = 0;
  }
  if (flagm == 1) {
    p = gsl_permutation_alloc(m->size1);
    gsl_linalg_complex_LU_decomp(mtmp, p, &signum);
  } else {
    if (itmp != argc-1) rb_raise(rb_eArgError, "signum not given");
    signum = NUM2DBL(argv[itmp]);
  }
  vz = Data_Make_Struct(cgsl_complex, gsl_complex, 0, free, z);
  *z = gsl_linalg_complex_LU_sgndet(mtmp, signum);
  if (flagm == 1) {
    gsl_matrix_complex_free(mtmp);
    gsl_permutation_free(p);
  }
  return vz;
}

static VALUE rb_gsl_linalg_complex_LU_invert(int argc, VALUE *argv, VALUE obj)
{
  gsl_matrix_complex *m = NULL, *mtmp = NULL, *inverse = NULL;
  gsl_permutation *p = NULL;
  int flagm = 0, signum, itmp;
  switch (TYPE(obj)) {
  case T_MODULE:
  case T_CLASS:
  case T_OBJECT:
    CHECK_MATRIX_COMPLEX(argv[0]);
    Data_Get_Struct(argv[0], gsl_matrix_complex, m);
    if (CLASS_OF(argv[0]) != cgsl_matrix_complex_LU) {
      mtmp = gsl_matrix_complex_alloc(m->size1, m->size2);
      gsl_matrix_complex_memcpy(mtmp, m);
      flagm = 1;
    } else {
      mtmp = m;
    }
    itmp = 1;
    break;
  default:
    Data_Get_Struct(obj, gsl_matrix_complex, m);
    if (CLASS_OF(obj) != cgsl_matrix_complex_LU) {
      mtmp = gsl_matrix_complex_alloc(m->size1, m->size2);
      gsl_matrix_complex_memcpy(mtmp, m);
      flagm = 1;
    } else {
      mtmp = m;
    }
    itmp = 0;
  }

  if (flagm == 1) {
    p = gsl_permutation_alloc(m->size1);
    gsl_linalg_complex_LU_decomp(mtmp, p, &signum);
  } else {
    Data_Get_Struct(argv[itmp], gsl_permutation, p);
  }
  inverse = gsl_matrix_complex_alloc(m->size1, m->size2);
  gsl_linalg_complex_LU_invert(mtmp, p, inverse);
  if (flagm == 1) {
    gsl_matrix_complex_free(mtmp);
    gsl_permutation_free(p);
  }
  return Data_Wrap_Struct(cgsl_matrix_complex, 0, gsl_matrix_complex_free, inverse);
}