matlisp-users

Re: [Matlisp-users] No Invert Functions?

[Michael A. K. <ma...@ll...>]

> I've just installed matlisp from the latest CVS at sourceforge. From 
> looking at the online help and reading some of the source, it looks to 
> me like none of the wrappers around matrix inversion functions (for 
> example, DGETRI and ZGETRI in LAPACK) have been implemented. Is there 
> any particular reason for this omission, or is it just lack of time? Or 
> am I totally missing something?

Joe,

  I don't see such a wrapper either.  However, matrix inversion is
currently performed with ZGESV by calling the MATLISP function M/ with
a single argument.  So if you just need the inverse of A use (M/ A).
Of course, if its ?GETRI specifically, a wrapper will be needed.


mike
**************************************************
Dr Michael A. Koerber    	Micro$oft Free Zone
MIT/Lincoln Laboratory	 	

Re: [Matlisp-users] No Invert Functions?

[Michael A. K. <ma...@ll...>]

> A related thing: I have recently implemented wrappers for
> dgetrs/zgetrs.  Are the maintainers interested in that change?  It
> involves changes in getrf.lisp, lapack.lisp, packages.lisp and
> system.dcl, as well as a new file getrs.lisp.  If yes, in which form
> do you want the contribution?

I'm just a user, however, I'd be interested in having this as part
of the MATLISP core.

mike
**************************************************
Dr Michael A. Koerber    	Micro$oft Free Zone
MIT/Lincoln Laboratory	 	

Re: [Matlisp-users] No Invert Functions?

[Nicolas N. <Nic...@iw...>]

Tunc Simsek <si...@ee...> writes:

> Hello; I can put the getrs stuff in the repository for you.
> It looks like you've got the code in the Matlisp style
> already so it should be straightforward.  I will need to know
> which functions are exported.  These functions should have
> doc strings.
> 
> I looked up getrs and it seems that it does basically m/ 
> but more explicitly.  That is, it keeps the residuals of the
> LU factorization.  Is there more to it?
> 
> Best, Tunc

The point is that computing the inverse is often not recommended
because of numerical stability.  Usually one uses the LU decomposition
directly which is done by the getrf/getrs-pair.  Of course, if you
only need to apply the inverse once you could be satisfied with m/.

Nicolas.


>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
Here are the patches (I used "cvs diff", is there a better way with
"cvs rdiff"?).  (I didn't change getrf! in contrast to what I said before.)
>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

~/CL-HOME/matlisp$ cvs diff system.dcl
Index: system.dcl
===================================================================
RCS file: /cvsroot/matlisp/matlisp/system.dcl,v
retrieving revision 1.19
diff -r1.19 system.dcl
183c183,184
< 		     "getrf"))
---
> 		     "getrf"
> 		     "getrs"))



~/CL-HOME/matlisp$ cvs diff packages.lisp
Index: packages.lisp
===================================================================
RCS file: /cvsroot/matlisp/matlisp/packages.lisp,v
retrieving revision 1.14
diff -r1.14 packages.lisp
118,119c118,119
<    "DGESV" "DGEEV" "DGETRF" "DGESVD"
<    "ZGESV" "ZGEEV" "ZGETRF" "ZGESVD" 
---
>    "DGESV" "DGEEV" "DGETRF" "DGETRS" "DGESVD"
>    "ZGESV" "ZGEEV" "ZGETRF" "ZGETRS" "ZGESVD" 
218a219,220
>      "GETRS!"
>      "GETRS"


And finally src/getrs.lisp:

;;; -*- Mode: lisp; Syntax: ansi-common-lisp; Package: :matlisp; Base: 10 -*-
;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;
;;; Copyright (c) 2000 The Regents of the University of California.
;;; All rights reserved. 
;;; 
;;; Permission is hereby granted, without written agreement and without
;;; license or royalty fees, to use, copy, modify, and distribute this
;;; software and its documentation for any purpose, provided that the
;;; above copyright notice and the following two paragraphs appear in all
;;; copies of this software.
;;; 
;;; IN NO EVENT SHALL THE UNIVERSITY OF CALIFORNIA BE LIABLE TO ANY PARTY
;;; FOR DIRECT, INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES
;;; ARISING OUT OF THE USE OF THIS SOFTWARE AND ITS DOCUMENTATION, EVEN IF
;;; THE UNIVERSITY OF CALIFORNIA HAS BEEN ADVISED OF THE POSSIBILITY OF
;;; SUCH DAMAGE.
;;;
;;; THE UNIVERSITY OF CALIFORNIA SPECIFICALLY DISCLAIMS ANY WARRANTIES,
;;; INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
;;; MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. THE SOFTWARE
;;; PROVIDED HEREUNDER IS ON AN "AS IS" BASIS, AND THE UNIVERSITY OF
;;; CALIFORNIA HAS NO OBLIGATION TO PROVIDE MAINTENANCE, SUPPORT, UPDATES,
;;; ENHANCEMENTS, OR MODIFICATIONS.
;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;
;;; Written by Nicolas Neuss (analogous to gesv.lisp by R. Toy)
;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;
;;; $Id:$
;;;
;;; $Log:$
;;; o Initial revision.
;;;
;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(in-package "MATLISP")

#+nil (use-package "BLAS")
#+nil (use-package "LAPACK")
#+nil (use-package "FORTRAN-FFI-ACCESSORS")

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;
;;; Generic function definitions
;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(defgeneric getrs! (a ipiv b &key trans)
  (:documentation
   "
  Syntax
  ======
  (GETRS! a ipiv b [:trans :N])

  Purpose
  =======
  Solves a system of linear equations
      A * X = B  or  A' * X = B
  with a general N-by-N matrix A using the LU factorization computed
  by GETRF.  A and IPIV are the results from GETRF, TRANS specifies
  the form of the system of equations: 
           = 'N':  A * X = B  (No transpose)
           = 'T':  A'* X = B  (Transpose)
           = 'C':  A'* X = B  (Conjugate transpose)

  Return Values
  =============
  [1] The NxM matrix X. (overwriting B)
  [4] INFO = T: successful
             i:  U(i,i) is exactly zero.  The LU factorization
                 used in the computation has been completed, 
                 but the factor U is exactly singular.
                 Solution could not be computed.
"))

(defgeneric getrs (a ipiv b &key trans)
  (:documentation
 "
  Sytnax
  ======
  (GETRS a ipiv b [:trans :N])

  Purpose
  =======
  Same as GETRS! except that B is not overwritten.
"))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;
;;; Method definitions
;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(defmethod getrs! :before ((a standard-matrix) ipiv (b standard-matrix) &key trans)
  (declare (ignore trans))
  (let ((n-a (nrows a))
	(m-a (ncols a))
	(n-b (nrows b)))
    (if (not (= n-a m-a n-b))
	(error "Dimensions of A,B given to GETRS! do not match"))
    ;;(check-type ipiv (simple-array (unsigned-byte 32) (*)))
    (if (< (length ipiv) n-a)
	(error "The argument IPIV given to GETRS! must have dimension >= N,
where NxN is the dimension of the argument A given to GETRS!"))))

(defmethod getrs! ((a real-matrix) ipiv (b real-matrix) &key trans)
  (let* ((n (nrows a))
	 (m (ncols b)))

    (declare (type fixnum n m)
	     (type (simple-array (unsigned-byte 32) (*)) ipiv))

    (multiple-value-bind (x info)
	(dgetrs (case trans
		  (:C "C")
		  (:T "T")
		  (t "N"))
		n
		m
		(store a)
		n
		ipiv
		(store b)
		n
		0)
	
	(values 
	 (make-instance 'real-matrix :nrows n :ncols m :store x)
	 (if (zerop info)
	     t
	   info)))))

(defmethod getrs! ((a complex-matrix) ipiv (b complex-matrix) &key trans)

  (let* ((n (nrows a))
	 (m (ncols b)))

    (declare (type fixnum n m)
	     (type (simple-array (unsigned-byte 32) (*)) ipiv))

    (multiple-value-bind (x info)
	(zgetrs (case trans
		  (:C "C")
		  (:T "T")
		  (t "N"))
		n
		m
		(store a)
		n
		ipiv
		(store b)
		n
		0)
	
	(values 
	 (make-instance 'complex-matrix :nrows n :ncols m :store x) 
	 (if (zerop info)
	     t
	   info)))))

(defmethod getrs! ((a standard-matrix) ipiv (b standard-matrix) &key trans)
  (let ((a (typecase a
	     (real-matrix (copy! a (make-complex-matrix-dim (nrows a) (ncols a))))
	     (complex-matrix a)
	     (t (error "argument A given to GETRS! is not a REAL-MATRIX or COMPLEX-MATRIX"))))
	(b (typecase b
	     (real-matrix (copy! b (make-complex-matrix-dim (nrows b) (ncols b))))
	     (complex-matrix b)
	     (t (error "argument B given to GETRS! is not a REAL-MATRIX or COMPLEX-MATRIX")))))
    (getrs! a ipiv b :trans trans)))

(defmethod getrs :before ((a standard-matrix) ipiv (b standard-matrix) &key trans)
  (declare (ignore trans))
  (if (not (= (nrows a) (ncols a) (nrows b) (length ipiv)))
      (error "dimensions of A,B,ipiv given to GETRS do not match")))

(defmethod getrs ((a real-matrix) ipiv (b real-matrix) &key trans)
  (getrs! a ipiv (copy b) :trans trans))

(defmethod getrs ((a complex-matrix) ipiv (b complex-matrix) &key trans)
  (getrs! a ipiv (copy b) :trans trans))

(defmethod getrs ((a standard-matrix) ipiv (b standard-matrix) &key trans)
  (let ((a (typecase a
	     (real-matrix (copy! a (make-complex-matrix-dim (nrows a) (ncols a))))
	     (complex-matrix (copy a))
	     (t (error "argument A given to GETRS is not a REAL-MATRIX or COMPLEX-MATRIX"))))
	(b (typecase b
	     (real-matrix (copy! b (make-complex-matrix-dim (nrows b) (ncols b))))
	     (complex-matrix (copy b))
	     (t (error "argument B given to GETRS is not a REAL-MATRIX or COMPLEX-MATRIX")))))
    (getrs! a ipiv b :trans trans)))

Re: [Matlisp-users] No Invert Functions?

[Nicolas N. <Nic...@iw...>]

I forgot:

cvs diff lapack.lisp
Index: lapack.lisp
===================================================================
RCS file: /cvsroot/matlisp/matlisp/src/lapack.lisp,v
retrieving revision 1.6
diff -r1.6 lapack.lisp
68,69c68,69
< "DGESV" "DGEEV" "DGETRF" "DGESVD"
< "ZGESV" "ZGEEV" "ZGETRF" "ZGESVD" ))
---
> "DGESV" "DGEEV" "DGETRF" "DGETRS" "DGESVD"
> "ZGESV" "ZGEEV" "ZGETRF" "ZGETRS" "ZGESVD" ))
321a322,386
> (def-fortran-routine dgetrs :void
> "
>    -- LAPACK routine (version 3.0) --
>       Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
>       Courant Institute, Argonne National Lab, and Rice University
>       March 31, 1993
> 
>    Purpose
>    =======
> 
>    DGETRS solves a system of linear equations
>       A * X = B  or  A' * X = B
>    with a general N-by-N matrix A using the LU factorization computed
>    by DGETRF.
> 
>    Arguments
>    =========
> 
>    TRANS   (input) CHARACTER*1
>            Specifies the form of the system of equations:
>            = 'N':  A * X = B  (No transpose)
>            = 'T':  A'* X = B  (Transpose)
>            = 'C':  A'* X = B  (Conjugate transpose = Transpose)
> 
>    N       (input) INTEGER
>            The order of the matrix A.  N >= 0.
> 
>    NRHS    (input) INTEGER
>            The number of right hand sides, i.e., the number of columns
>            of the matrix B.  NRHS >= 0.
> 
>    A       (input) DOUBLE PRECISION array, dimension (LDA,N)
>            The factors L and U from the factorization A = P*L*U
>            as computed by DGETRF.
> 
>    LDA     (input) INTEGER
>            The leading dimension of the array A.  LDA >= max(1,N).
> 
>    IPIV    (input) INTEGER array, dimension (N)
>            The pivot indices from DGETRF; for 1<=i<=N, row i of the
>            matrix was interchanged with row IPIV(i).
> 
>    B       (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
>            On entry, the right hand side matrix B.
>            On exit, the solution matrix X.
> 
>    LDB     (input) INTEGER
>            The leading dimension of the array B.  LDB >= max(1,N).
> 
>    INFO    (output) INTEGER
>            = 0:  successful exit
>            < 0:  if INFO = -i, the i-th argument had an illegal value
> 
> "
>   (trans :string :input)
>   (n :integer :input)
>   (nhrs :integer :input)
>   (a (* :double-float) :input)
>   (lda :integer :input)
>   (ipiv (* :integer) :input)
>   (b (* :double-float) :input-output)
>   (ldb :integer :input)
>   (info :integer :output)
> )
> 
694a760,823
>   (info :integer :output)
> )
> 
> (def-fortran-routine zgetrs :void
> "
>    -- LAPACK routine (version 3.0) --
>       Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
>       Courant Institute, Argonne National Lab, and Rice University
>       September 30, 1994
>  
>    Purpose
>    =======
>  
>    ZGETRS solves a system of linear equations
>       A * X = B,  A**T * X = B,  or  A**H * X = B
>    with a general N-by-N matrix A using the LU factorization computed
>    by ZGETRF.
>  
>    Arguments
>    =========
>  
>    TRANS   (input) CHARACTER*1
>            Specifies the form of the system of equations:
>            = 'N':  A * X = B     (No transpose)
>            = 'T':  A**T * X = B  (Transpose)
>            = 'C':  A**H * X = B  (Conjugate transpose)
>  
>    N       (input) INTEGER
>            The order of the matrix A.  N >= 0.
>  
>    NRHS    (input) INTEGER
>            The number of right hand sides, i.e., the number of columns
>            of the matrix B.  NRHS >= 0.
>  
>    A       (input) COMPLEX*16 array, dimension (LDA,N)
>            The factors L and U from the factorization A = P*L*U
>            as computed by ZGETRF.
>  
>    LDA     (input) INTEGER
>            The leading dimension of the array A.  LDA >= max(1,N).
>  
>    IPIV    (input) INTEGER array, dimension (N)
>            The pivot indices from ZGETRF; for 1<=i<=N, row i of the
>            matrix was interchanged with row IPIV(i).
>  
>    B       (input/output) COMPLEX*16 array, dimension (LDB,NRHS)
>            On entry, the right hand side matrix B.
>            On exit, the solution matrix X.
>  
>    LDB     (input) INTEGER
>            The leading dimension of the array B.  LDB >= max(1,N).
>  
>    INFO    (output) INTEGER
>            = 0:  successful exit
>            < 0:  if INFO = -i, the i-th argument had an illegal value
> "
>   (trans :string :input)
>   (n :integer :input)
>   (nhrs :integer :input)
>   (a (* :complex-double-float) :input)
>   (lda :integer :input)
>   (ipiv (* :integer) :input)
>   (ldb :integer :input)
>   (b (* :complex-double-float) :input-output)

Re: [Matlisp-users] No Invert Functions?

[Michael A. K. <ma...@ll...>]

Nicolas,

1.  (M/ A) eventually calls (GESV A I), where I is the identity
    matrix, which intern calls the FORTRAN ?GESV which performs the LU
    decomposition on the system of equations [A I] to form the inverse
    of A.

2.  Your statement concerning inverses is correct...and it is correctly 
    implimented in MATLISP.

> The point is that computing the inverse is often not recommended
> because of numerical stability.  Usually one uses the LU decomposition
> directly which is done by the getrf/getrs-pair.  Of course, if you
> only need to apply the inverse once you could be satisfied with m/.
> 


mike
**************************************************
Dr Michael A. Koerber    	Micro$oft Free Zone
MIT/Lincoln Laboratory	 	

Re: [Matlisp-users] No Invert Functions?

[Nicolas N. <Nic...@iw...>]

"Michael A. Koerber" <ma...@ll...> writes:

> Nicolas,
> 
> 1.  (M/ A) eventually calls (GESV A I), where I is the identity
>     matrix, which intern calls the FORTRAN ?GESV which performs the LU
>     decomposition on the system of equations [A I] to form the inverse
>     of A.
> 
> 2.  Your statement concerning inverses is correct...and it is correctly 
>     implimented in MATLISP.
> 
> > The point is that computing the inverse is often not recommended
> > because of numerical stability.  Usually one uses the LU decomposition
> > directly which is done by the getrf/getrs-pair.  Of course, if you
> > only need to apply the inverse once you could be satisfied with m/.
> > 


The problem with using only GESV and M/ is that it computes an LU
decomposition (which is the major O(N^3) part of the work).  To use it
again (e.g. in a time-stepping scheme) you need GETRS.

Nicolas.

Re: [Matlisp-users] No Invert Functions?

[Nicolas N. <Nic...@iw...>]

Nicolas Neuss <Nic...@iw...> writes:

> The problem with using only GESV and M/ is that it computes an LU
> decomposition (which is the major O(N^3) part of the work).

This is written badly.  The problem is that this decomposition cannot
be used again without GETRS.

>  To use it again (e.g. in a time-stepping scheme) you need GETRS.
> 
> Nicolas.

Re: [Matlisp-users] No Invert Functions?

[Michael A. K. <ma...@ll...>]

> The problem with using only GESV and M/ is that it computes an LU
> decomposition (which is the major O(N^3) part of the work).  To use it
> again (e.g. in a time-stepping scheme) you need GETRS.
> 
Nicolas,

1.  I have no argument with this...I do have argument with implying that
    MATLISP does _not_ now use LU decomposition to compute inverses.

mike

Re: [Matlisp-users] pseudoinverse

[Michael A. K. <ma...@ll...>]

Nicolas,

> I don't know about the other members of this mailing list, but what
> concerns me: I would like to see discussions about new add-ons on the
> matlisp-user list.  The traffic is not large here and it would be nice

My initial post, to mat...@li... is in compliance
with your request.  The only thing more I could have done was to post
my thought process...not something that most want to see :-).

So...disucssion is now open...do Matlisp user want a PINV command added?


**************************************************
Dr Michael A. Koerber    	Micro$oft Free Zone
MIT/Lincoln Laboratory	 	

Re: [Matlisp-users] pseudoinverse

[Tunc S. <si...@ee...>]

Hi Michael,

Thanks so much for PINV.  It solves a basic problem, and
I think it will make a useful
addition to Matlisp.  Just one more thing to think about.
Can you say something about the difference between
your approach (using SVD) and the LAPACK routines called
DGELS and DGELSD. Numerical stability?

Your question on MACH pinpoints a fundamental issue in Matlisp.
Almost all (if not all) Matlisp routines use the underlying
fortran implementation to do the work.  Unless anything has
changed, the MACH is the one that fortran would use by default.
If this is not sufficient, or you need explicit access to this
number, we need to do some extra work (most likely in the
config file for matlisp). 

Best, Tunc



"Michael A. Koerber" wrote:
> 
> Nicolas,
> 
> > I don't know about the other members of this mailing list, but what
> > concerns me: I would like to see discussions about new add-ons on the
> > matlisp-user list.  The traffic is not large here and it would be nice
> 
> My initial post, to mat...@li... is in compliance
> with your request.  The only thing more I could have done was to post
> my thought process...not something that most want to see :-).
> 
> So...disucssion is now open...do Matlisp user want a PINV command added?
> 
> **************************************************
> Dr Michael A. Koerber           Micro$oft Free Zone
> MIT/Lincoln Laboratory
> 
> -------------------------------------------------------
> This sf.net email is sponsored by:ThinkGeek
> Welcome to geek heaven.
> http://thinkgeek.com/sf
> _______________________________________________
> Matlisp-users mailing list
> Mat...@li...
> https://lists.sourceforge.net/lists/listinfo/matlisp-users

Re: [Matlisp-users] No Invert Functions?

[Joseph D. <jd...@uc...>]

Michael A. Koerber wrote:
>>I've just installed matlisp from the latest CVS at sourceforge. From 
>>looking at the online help and reading some of the source, it looks to 
>>me like none of the wrappers around matrix inversion functions (for 
>>example, DGETRI and ZGETRI in LAPACK) have been implemented. Is there 
>>any particular reason for this omission, or is it just lack of time? Or 
>>am I totally missing something?
> 
> 
> Joe,
> 
>   I don't see such a wrapper either.  However, matrix inversion is
> currently performed with ZGESV by calling the MATLISP function M/ with
> a single argument.  So if you just need the inverse of A use (M/ A).
> Of course, if its ?GETRI specifically, a wrapper will be needed.
> 

Thanks Michael and Tunc,

M/ should be sufficient for me. However, there appears to be a bug/typo 
in the documentation:


* (help m/)

In package MATLISP, the function M/
===================================

   Syntax
   ======
   (M/ a [b])

   Purpose
   =======

   Equivalent to inv(A) * B, but far more efficient and accurate.
   If A is not given, the inverse of B is returned.

   A is an NxN square matrix and B is NxM.

   B may be a scalar in which case M/ is equivalent to M./

   See M/!, M./, M./!, GESV


It seems that the line "If A is not given, the inverse of B is 
returned." should read "If B is not given, the inverse of A is returned."


Joe

Re: [Matlisp-users] No Invert Functions?

[Raymond T. <to...@rt...>]

>>>>> "Joseph" == Joseph Dale <jd...@uc...> writes:

    Joseph> Michael A. Koerber wrote:
    >>> I've just installed matlisp from the latest CVS at
    >>> sourceforge. From looking at the online help and reading some of
    >>> the source, it looks to me like none of the wrappers around matrix
    >>> inversion functions (for example, DGETRI and ZGETRI in LAPACK) have
    >>> been implemented. Is there any particular reason for this omission,
    >>> or is it just lack of time? Or am I totally missing something?

    >> Joe,

    >> I don't see such a wrapper either.  However, matrix inversion is

    >> currently performed with ZGESV by calling the MATLISP function M/ with
    >> a single argument.  So if you just need the inverse of A use (M/ A).
    >> Of course, if its ?GETRI specifically, a wrapper will be needed.
    >> 


    Joseph> Thanks Michael and Tunc,

    Joseph> M/ should be sufficient for me. However, there appears to be a
    Joseph> bug/typo in the documentation:

Thanks.

I'll fix this.

Ray