Class: NMatrixLU

Inherits:
Object
  • Object
show all
Defined in:
ext/narray/na_linalg.c

Instance Method Summary collapse

Constructor Details

#initialize(lu, piv) ⇒ Object

:nodoc:



543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
 
# File 'ext/narray/na_linalg.c', line 543

static VALUE
na_lu_init(VALUE self, VALUE lu, VALUE piv)
{
  int i;
  struct NARRAY *l, *p;

  if (CLASS_OF(lu)!=cNMatrix)
    rb_raise(rb_eTypeError,"LU should be NMatrix");
  if (CLASS_OF(piv)!=cNVector)
    rb_raise(rb_eTypeError,"pivot should be NVector");

  GetNArray(lu,l);
  GetNArray(piv,p);

  if (p->type != NA_LINT)
    rb_raise(rb_eRuntimeError,"pivot type must be Integer");

  if (l->rank != p->rank+1)
    rb_raise(rb_eRuntimeError,"array dimension mismatch %i!=%i+1",
	     l->rank, p->rank);

  if (l->shape[0] != l->shape[1])
    rb_raise(rb_eRuntimeError,"LU matrix (%zd,%zd) is not square",
	     l->shape[0], l->shape[1]);

  for (i=1; i<l->rank; ++i)
    if (l->shape[i] != p->shape[i-1])
      rb_raise(rb_eRuntimeError,"array size mismatch %zd!=%zd at %i",
	       l->shape[i], p->shape[i-1], i);

  rb_ivar_set(self, id_lu, lu);
  rb_ivar_set(self, id_pivot, piv);
  return Qnil;
}

Instance Method Details

#solve(arg) ⇒ Object

Solve with the result of LU factorization. arg should be NMatrix or NVector instance. Returns an instance of same class with arg.



487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
 
# File 'ext/narray/na_linalg.c', line 487

static VALUE
na_lu_solve(VALUE self, volatile VALUE other)
{
  int  ndim;
  na_shape_t n, *shape;
  struct NARRAY *a1, *a2, *l, *p;
  VALUE pv, obj, klass;
  volatile VALUE lu;

  klass = CLASS_OF(other);
  if (klass==cNVector)
    other = na_newdim_ref(1,(VALUE*)na_funcset[NA_ROBJ].zero,other);
  else if (klass!=cNMatrix)
    rb_raise(rb_eTypeError,"neither NMatrix or NVector");

  lu = rb_ivar_get(self, id_lu);
  pv = rb_ivar_get(self, id_pivot);

  GetNArray(lu,l);

  other = na_upcast_object(other,l->type);
  GetNArray(other,a1);

  lu = na_upcast_type(lu,a1->type);
  GetNArray(lu,l);
  GetNArray(pv,p);

  n = l->shape[0];
  if (n != a1->shape[1])
    rb_raise(rb_eTypeError,"size mismatch (%zd!=%zd)",n,a1->shape[1]);

  ndim  = NA_MAX(l->rank, a1->rank);
  shape = ALLOCA_N(na_shape_t, ndim);

  shape[0] = a1->shape[0];
  na_shape_max2( ndim-1, shape+1, a1->rank-1, a1->shape+1,
		 l->rank-1, l->shape+1 );
  obj = na_make_object( a1->type, ndim, shape, klass );

  GetNArray(obj,a2);

  na_exec_linalg( a2, a1, p, 2, 2, 1, na_lu_pivot_func );
  na_exec_linalg( a2, a2, l, 2, 2, 2, na_lu_solve_func );

  if (klass==cNVector) {
    shape = ALLOC_N(na_shape_t, ndim-1);
    memcpy(shape,a2->shape+1,sizeof(na_shape_t)*(ndim-1));
    xfree(a2->shape);
    a2->shape = shape;
    --(a2->rank);
  }
  return obj;
}