[Matlisp-users] Cholesky decomp. and sol.

[Knut G. <knu...@nt...>]

Hello,

I want to be able to use the Cholesky method in matlisp, i.e. access the
functions dpotrf and dpotrs from lapack. I wrote the files needed, basically
copying from getrf.lisp and getrs.lisp in matlisp/src. It compiles and
everything looks OK. However, running 

(let* ((a (matlisp:make-real-matrix-dim 3 3))(ipiv (make-array 3 :element-type
'(unsigned-byte 32)))(b (matlisp:make-real-matrix-dim 3 1)))(loop for i from 0
to 2 do (loop for j from 0 to 2 do (setf (matlisp:matrix-ref a i j)
10.0d0)))(setf (matlisp:matrix-ref a 0 0) 25.0d0)(setf (matlisp:matrix-ref a 2
2) 62.0d0)(setf (matlisp:matrix-ref a 0 1) -5.0d0)(setf (matlisp:matrix-ref a 1
0) -5.0d0)(setf (matlisp:matrix-ref a 1 1) 17.0d0)(setf (matlisp:matrix-ref b 0
0) 55.0d0)(setf (matlisp:matrix-ref b 1 0) -19.0d0)(setf (matlisp:matrix-ref b 2
0) 114.0d0)(matlisp:getrf! a ipiv )(matlisp:getrs! a ipiv b :trans :N))

gives the correct solution [1,-2,2], whilst 

(let* ((a (matlisp:make-real-matrix-dim 3 3))(b (matlisp:make-real-matrix-dim 3
1 )))(loop for i from 0 to 2 do (loop for j from 0 to 2 do (setf
(matlisp:matrix-ref a i j) 10.0d0)))(setf (matlisp:matrix-ref a 0 0)
25.0d0)(setf (matlisp:matrix-ref a 2 2) 62.0d0)(setf (matlisp:matrix-ref a 0 1)
-5.0d0)(setf (matlisp:matrix-ref a 1 0) -5.0d0)(setf (matlisp:matrix-ref a 1 1)
17.0d0)(setf (matlisp:matrix-ref b 0 0) 55.0d0)(setf (matlisp:matrix-ref b 1 0)
-19.0d0)(setf (matlisp:matrix-ref b 2 0) 114.0d0)(matlisp:potrf! a :uplo
:U)(print a)(print b)(matlisp:potrs! a b :uplo :U))

gives an error message:

#<MATLISP:REAL-MATRIX  3 x 3
    5.0000     -1.0000      2.0000     
   -5.0000      4.0000      3.0000     
    10.000      10.000      7.0000     > 
#<MATLISP:REAL-MATRIX  3 x 1
    55.000     
   -19.000     
    114.00     > 
debugger invoked on a SIMPLE-TYPE-ERROR:
  The value of FORTRAN-FFI-ACCESSORS::VEC is 3, which is not of type
FORTRAN-FFI-ACCESSORS::MATLISP-SPECIALIZED-ARRAY.

Type HELP for debugger help, or (SB-EXT:QUIT) to exit from SBCL.

restarts (invokable by number or by possibly-abbreviated name):
  0: [STORE-VALUE] Supply a new value for FORTRAN-FFI-ACCESSORS::VEC.
  1: [ABORT      ] Exit debugger, returning to top level.

(SB-KERNEL:CHECK-TYPE-ERROR
 FORTRAN-FFI-ACCESSORS::VEC
 3
 FORTRAN-FFI-ACCESSORS::MATLISP-SPECIALIZED-ARRAY
 NIL)
0] 
WARNING: Starting a select without a timeout while interrupts are disabled.

I threw in the (print )s just to verify that both matrices are of type
real-matrix. Why does this not give the correct answer and wat do I have to do?

Thanks,

Knut Gjerden