matlisp-users — General discussion list.

[Matlisp-users] allegro heap size limit?

[Jefferson P. <jp...@cs...>]

This is not strictly a MATLISP question, but I though people on this
list might have some experience with very large Lisp heaps.

I'm using Allegro 6.0, on a P4 machine running RedHat with 1.5GB of RAM 
and 3GB of swap.  I'm unable to get a heap-size larger than about 500MB. 
  As soon as an allocation attempts to grow the heap above that level, I 
get something like this:

---
Error: An allocation request for 4016 bytes caused a need for 171442176 
more bytes of heap.  The operating system will not make the space
available because of a lack of swap space or some other operating system
imposed limit or memory mapping collision.
   [condition type: STORAGE-CONDITION]
---

When I got this error the lisp process size was about 360MB.  Allocating 
171MB would have put it at 530MB.

A C program on the same machine has no problem allocating blocks of 1GB
or more, which makes me wonder if it's really the OS generating the error...

J.

Re: [Matlisp-users] CMUCL 18d failing to compile...

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

>>>>> "Jeremy" == Jeremy Shute <sh...@op...> writes:

    Jeremy> Hmmm, I see that in neither 1.0A or 1.0B's install files.  I tried
    Jeremy> building 1.0A tonight, and I'm getting something similar:

    Jeremy> ("DEFUN MAKE-PATHNAME" 268433206 8)[:OPTIONAL]
    Jeremy> Source: Error finding source:
    Jeremy> Error in function DEBUG::GET-FILE-TOP-LEVEL-FORM:  Source file no longer
    Jeremy> exists:  target:code/pathname.lisp.

    Jeremy> ...did you build from CVS?

Yes, please build from CVS.  The releases are hopelessly out-of-date.
My apologies for that.

Ray

Re: [Matlisp-users] CMUCL 18d failing to compile...

[Jeremy S. <sh...@op...>]

Hmmm, I see that in neither 1.0A or 1.0B's install files.  I tried
building 1.0A tonight, and I'm getting something similar:

("DEFUN MAKE-PATHNAME" 268433206 8)[:OPTIONAL]
Source: Error finding source:
Error in function DEBUG::GET-FILE-TOP-LEVEL-FORM:  Source file no longer
exists:  target:code/pathname.lisp.

...did you build from CVS?

Jeremy



On Wed, 2003-03-12 at 07:48, Michael A. Koerber wrote:
> In the INSTALL file read that the sequence of events is something like
> 
>   ./configure --with-lisp=cmucl
>   make
> 
> 
> > 4. Configure the system.
> >    o If you're using Allegro CL, specify --with-lisp=acl.
> >    o If you're using CMU CL, specify --with-lisp=cmucl
> >    o Give the name of the Lisp executable via --with-lisp-exec=<name>
> >    o Run configure:
> > 
> >    ./configure --with-lisp=<lisp> --with-lisp-exec=<exec> --prefix=`pwd`
> > 
> >    If you are not using g77 to compile, you may want to specify the
> >    Fortran flags.  For example, to build with Sun's Fortran compiler
> >    this can be used (from bash):
> > 
> >    F77=f77 FFLAGS='-g -O -KPIC' ./configure ...
> > 
> > 5. Build the Fortran and Lisp files
> > 
> > 	make
> 
> 
> mike
> 

Re: [Matlisp-users] CMUCL 18d failing to compile...

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

In the INSTALL file read that the sequence of events is something like

  ./configure --with-lisp=cmucl
  make


> 4. Configure the system.
>    o If you're using Allegro CL, specify --with-lisp=acl.
>    o If you're using CMU CL, specify --with-lisp=cmucl
>    o Give the name of the Lisp executable via --with-lisp-exec=<name>
>    o Run configure:
> 
>    ./configure --with-lisp=<lisp> --with-lisp-exec=<exec> --prefix=`pwd`
> 
>    If you are not using g77 to compile, you may want to specify the
>    Fortran flags.  For example, to build with Sun's Fortran compiler
>    this can be used (from bash):
> 
>    F77=f77 FFLAGS='-g -O -KPIC' ./configure ...
> 
> 5. Build the Fortran and Lisp files
> 
> 	make


mike

[Matlisp-users] CMUCL 18d failing to compile...

[Jeremy S. <sh...@op...>]

Hi:

I'm trying to compile Matlisp using the RPM install of CMUCL 18d; seems
like I'm not getting the install correct...  I'm very new to LISP
(though I'm quite familiar with Scheme), adding it to the repertoire
specifically for this library.  (Sick of writing things in Octave and
CLAPACK.)

Unpacked and ran 'make cmu', didn't fiddle with 'configure' as the
installation notes made no mention of it.

Following is a typescript of what's going on...

Script started on Tue Mar 11 21:23:52 2003
[shutej@localhost matlisp-1.0b]$ make cmu
/usr/local/bin/lisp -eval '(progn (load "start.lisp"))'
; Loading #p"/home/shutej/tmp/matlisp-1.0b/start.lisp".
;; Loading #p"/home/shutej/tmp/matlisp-1.0b/config.lisp".
;; Loading #p"/home/shutej/tmp/matlisp-1.0b/system.dcl".
;;; Loading #p"/home/shutej/tmp/matlisp-1.0b/packages.lisp".
;;; Loading #p"/home/shutej/tmp/matlisp-1.0b/packages.lisp".
;;; Loading #p"/home/shutej/tmp/matlisp-1.0b/packages.lisp".
;;; Loading #p"/home/shutej/tmp/matlisp-1.0b/packages.lisp".
;;; Loading #p"/home/shutej/tmp/matlisp-1.0b/packages.lisp".
;;; Loading #p"/home/shutej/tmp/matlisp-1.0b/defsystem.lisp".

Warning: These variables are undefined:
  *LIBRARY* *MODULE-FILES*
Warning:  Old-style IN-PACKAGE.

In: LAMBDA (FILENAME &REST ARGS &KEY OUTPUT-FILE ...)
  #'(LAMBDA (FILENAME &REST ARGS &KEY OUTPUT-FILE ...)
      (DECLARE (IGNORE ARGS))
      (BLOCK C-COMPILE-FILE (RUN-UNIX-PROGRAM *C-COMPILER* `#)))
Note: Variable ERROR-FILE defined but never used.

Type-error in KERNEL::OBJECT-NOT-TYPE-ERROR-HANDLER:
   #() is not of type (OR CONS
                          BASE-STRING
                          (MEMBER NIL :UNSPECIFIC :WILD)
                          COMMON-LISP::PATTERN)

Restarts:
  0: [CONTINUE] Return NIL from load of "matlisp:system.dcl".
  1:            Return NIL from load of "start.lisp".
  2: [ABORT   ] Skip remaining initializations.

Debug  (type H for help)

("DEFUN MAKE-PATHNAME" 268433174 8)[:OPTIONAL]
Source: Error finding source: 
Error in function DEBUG::GET-FILE-TOP-LEVEL-FORM:  Source file no longer
exists:
  target:code/pathname.lisp.
0] (quit)
[shutej@localhost matlisp-1.0b]$ exit

Script done on Tue Mar 11 21:24:01 2003

...CMUCL came from the RPMs linked to from the CMUCL website (with the
Extras installed as well.)  Redhat Linux 8.0 running on an Intel P4.

Thanks,

Jeremy

[Matlisp-users] 0x0-matrices

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

Hello,

at the moment matlisp can't handle 0x0-, nx0- or 0xn-matrices.  This is
annoying because these may arise in some border situations (of course,
without elements being referenced).

Below is a patch which should make this possible (it corrects also a
remaining :type 'fixnum).  Sourceforge CVS was partially down, so I diffed
against a local file for print.lisp (which should be same as CVS, though).

Yours, Nicolas.


neuss@ortler:~/CL-HOME/matlisp/src$ cvs diff -u matrix.lisp
Index: matrix.lisp
===================================================================
RCS file: /cvsroot/matlisp/matlisp/src/matrix.lisp,v
retrieving revision 1.11
diff -u -r1.11 matrix.lisp
--- matrix.lisp	19 Feb 2003 21:59:52 -0000	1.11
+++ matrix.lisp	28 Feb 2003 15:20:19 -0000
@@ -256,7 +256,7 @@
     :initarg :store-size
     :initform 0
     :accessor store-size
-    :type 'fixnum
+    :type fixnum
     :documentation "Total number of elements needed to store the matrix.  (Usually
 the same as nels, but not necessarily so!")
    (store
@@ -665,7 +665,7 @@
        (let ((arg (first args)))
 	 (typecase arg
 	   (integer
-	    (assert (plusp arg) nil
+	    (assert (not (minusp arg)) nil
 		    "matrix dimension must be positive, not ~A" arg)
 	    (make-real-matrix-dim arg arg))
 	   (sequence
@@ -676,8 +676,8 @@
       (2
        (destructuring-bind (n m)
 	   args
-	 (assert (and (typep n '(integer 1))
-		      (typep n '(integer 1)))
+	 (assert (and (typep n '(integer 0))
+		      (typep n '(integer 0)))
 		 nil
 		 "cannot make a ~A x ~A matrix" n m)
 	 (make-real-matrix-dim n m)))
@@ -850,8 +850,8 @@
        (let ((arg (first args)))
 	 (typecase arg
 	   (integer
-	    (assert (plusp arg) nil
-		    "matrix dimension must be positive, not ~A" arg)
+	    (assert (not (minusp arg)) nil
+		    "matrix dimension must be non-negative, not ~A" arg)
 	    (make-complex-matrix-dim arg arg))
 	   (sequence
 	    (make-complex-matrix-sequence arg))
@@ -861,8 +861,8 @@
       (2
        (destructuring-bind (n m)
 	   args
-	 (assert (and (typep n '(integer 1))
-		      (typep n '(integer 1)))
+	 (assert (and (typep n '(integer 0))
+		      (typep n '(integer 0)))
 		 nil
 		 "cannot make a ~A x ~A matrix" n m)
 	 (make-complex-matrix-dim n m)))

neuss@ortler:~/CL-HOME/matlisp/src$ diff -u print.lisp print.lisp.~1.7~
--- print.lisp	Fri Feb 28 16:30:46 2003
+++ print.lisp.~1.7~	Fri Feb 28 16:32:20 2003
@@ -154,35 +154,36 @@
 		  (dotimes (i *matrix-indent*)
 		    (format stream " "))
 		  (dotimes (j max-m)
+		    (declare (type fixnum j))
 		    (print-element matrix 
 				   (matrix-ref matrix i j)
 				   stream)
 		    (format stream " "))
-		  (when (plusp number-of-cols)
-		    (if (< max-m (1- number-of-cols))
+		  (if (< max-m (1- number-of-cols))
+		      (progn
+			(format stream "... ")
+			(print-element matrix 
+				       (matrix-ref matrix i (1- number-of-cols))
+				       stream)
+			(format stream " "))
+		    (if (< max-m number-of-cols)
 			(progn
-			  (format stream "... ")
 			  (print-element matrix 
 					 (matrix-ref matrix i (1- number-of-cols))
 					 stream)
-			  (format stream " "))
-			(if (< max-m number-of-cols)
-			    (progn
-			      (print-element matrix 
-					     (matrix-ref matrix i (1- number-of-cols))
-					     stream)
-			      (format stream " ")))))))
+			  (format stream " "))))))
 	   
 	   (dotimes (i max-n)
+	     (declare (type fixnum i))
 	     (print-row i))
+	   
+	   (if (< max-n (1- number-of-rows))
+	       (progn
+		 (format stream "~%     :")
+		 (print-row (1- number-of-rows)))
+	     (if (< max-n number-of-rows)
+		 (print-row (1- number-of-rows))))))))
 
-	   (when (plusp number-of-rows)
-	     (if (< max-n (1- number-of-rows))
-		 (progn
-		   (format stream "~%     :")
-		   (print-row (1- number-of-rows)))
-		 (if (< max-n number-of-rows)
-		     (print-row (1- number-of-rows)))))))))
 
 (defmethod print-object ((matrix standard-matrix) stream)
   (print-unreadable-object (matrix stream :type t :identity (not *print-matrix*))

Re: [Matlisp-users] FFT/IFFT Seriously Broke for Matrices

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

>>>>> "Mike" == mak  <ma...@ll...> writes:

    Mike> Ray,
    Mike>  tnx...I've made some minor (cosmetic) changes to the fft.lisp source.  
    Mike> If you care for them, they are below

I assume this should be used in place of your previous patch?

It seems to be ready to go so I'll check it in soon.

Ray

Re: [Matlisp-users] FFT/IFFT Seriously Broke for Matrices

[<ma...@ll...>]

Ray,

 tnx...I've made some minor (cosmetic) changes to the fft.lisp source.  
If you care for them, they are below

tnx

mike

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

*** fft.lisp.orig	Thu Oct 25 17:48:45 2001
--- fft.lisp	Tue Feb  4 11:45:16 2003
***************
*** 217,238 ****
  		       ((col-vector-p x)
  			(make-complex-matrix-dim n 1))
  		       (t (make-complex-matrix-dim n (ncols x))))))
      (if (row-or-col-vector-p x)
  	(progn
  	  (copy! x result)
! 	  (zfftf n (store result) wsave))
  
!       (dotimes (j (ncols x))
! 	(declare (type fixnum j))
! 	 (dotimes (i (nrows x))
! 	   (declare (type fixnum i))
! 	   (setf (matrix-ref result i j) (matrix-ref x i j)))
! 	 (with-vector-data-addresses ((addr-result (store result))
! 				      (addr-wsave wsave))
! 	    (incf-sap :complex-double-float addr-result (* j n))
! 	    (dfftpack::fortran-zfftf n addr-result addr-wsave))))
  
!       result))
  
  
  #+:allegro
--- 217,242 ----
  		       ((col-vector-p x)
  			(make-complex-matrix-dim n 1))
  		       (t (make-complex-matrix-dim n (ncols x))))))
+ 
      (if (row-or-col-vector-p x)
  	(progn
  	  (copy! x result)
! 	  (zfftf n (store result) wsave)
! 	  result)
  
!       ;; Do FFT by column...first we copy the array
!       (progn
! 	;; copy X to a working array, RESULT
! 	(dotimes (j-copy (ncols x))
! 	  (declare (type fixnum j-copy))
! 	  
! 	  ;; Note: we may be asked to truncate if N < (NROWS X)
! 	  (dotimes (i-copy (min n (nrows x))) 
! 	    (declare (type fixnum i-copy))
! 	    (setf (matrix-ref result i-copy j-copy) (matrix-ref x i-copy j-copy))))
  
! 	;; Then do FFT by column
! 	(fft! result)))))
  
  
  #+:allegro
***************
*** 286,313 ****
  	(progn
  	  (copy! x result)
  	  (zfftb n (store result) wsave)
! 	  (let ((scale-factor (/ (float n 1d0))))
! 	    (dotimes (k n)
! 	      (setf (matrix-ref result k)
! 		    (* scale-factor (matrix-ref result k))))))
! 
! 	(let ((scale-factor (/ (float n 1d0))))
! 	  (dotimes (j (ncols x))
! 	    (declare (type fixnum j))
! 	    (dotimes (i (nrows x))
! 	      (declare (type fixnum i))
! 	      (setf (matrix-ref result i j) (matrix-ref x i j)))
! 	    (with-vector-data-addresses ((addr-result (store result))
! 					 (addr-wsave wsave))
! 	      (incf-sap :complex-double-float addr-result (* j n))
! 	      (dfftpack::fortran-zfftb n addr-result addr-wsave))
! 	    ;; Scale the result
! 	    (dotimes (i (nrows x))
! 	      (declare (type fixnum i))
! 	      (setf (matrix-ref result i j) (* scale-factor (matrix-ref x i j)))))))
  
!     result))
  
  
  #+:allegro
  (defmethod ifft ((x standard-matrix) &optional n)
--- 290,310 ----
  	(progn
  	  (copy! x result)
  	  (zfftb n (store result) wsave)
! 	  (scal! (/ (float n 1d0)) result))
  
!       ;; Perform the IFFT by columns...first we copy the array
!       (progn
! 	;; Copy X to the working Array RESULT
! 	(dotimes (j-copy (ncols x))
! 	  (declare (type fixnum j-copy))
! 
! 	  ;; Note: we may be asked to truncate if N < (NROWS X)
! 	  (dotimes (i-copy (min n (nrows x)))
! 	    (declare (type fixnum i-copy))
! 	    (setf (matrix-ref result i-copy j-copy) (matrix-ref x i-copy j-copy))))
  
+ 	;; The we do the IFFT
+ 	(ifft! result)))))
  
  #+:allegro
  (defmethod ifft ((x standard-matrix) &optional n)
***************
*** 348,350 ****
--- 345,401 ----
  
        result))
  
+ ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+ ;; Modification for destructive FFT's.  The purpose is to eliminate
+ ;; cons'ing for certain application requiring numerous calls to these
+ ;; routines.
+ 
+ (defgeneric fft! (x)
+   (:documentation "See FFT but note that the optional N is NOT permitted.
+   Performs in place FFT modifying X.  This will only work for a complex matrix."))
+ 
+ (defgeneric ifft! (x)
+   (:documentation "See IFFT but note that the optional N is NOT permitted.
+   Performs in place IFFT modifying X.  This will only work for a complex matrix."))
+ 
+ #+:cmu  
+ (defmethod fft! ((x complex-matrix))
+   (let* ((n (if (row-or-col-vector-p x)
+ 		(max (nrows x) (ncols x))
+ 	      (nrows x)))
+ 	 (wsave (lookup-wsave-entry n)))
+ 
+     (if (row-or-col-vector-p x)
+ 	(progn
+ 	  (zfftf n (store x) wsave)
+ 	  x)
+ 
+       (dotimes (j (ncols x) x)
+ 	(declare (type fixnum j))
+ 	(with-vector-data-addresses ((addr-x (store x))
+ 				     (addr-wsave wsave))
+ 	    (incf-sap :complex-double-float addr-x (* j n))
+ 	    (dfftpack::fortran-zfftf n addr-x addr-wsave))))))
+ 
+ 
+ #+:cmu
+ (defmethod ifft! ((x complex-matrix))
+   (let* ((n (if (row-or-col-vector-p x)
+ 		      (max (nrows x) (ncols x))
+ 		      (nrows x)))
+ 	 (wsave (lookup-wsave-entry n)))
+ 
+     (if (row-or-col-vector-p x)
+ 	(progn
+ 	  (zfftb n (store x) wsave)
+ 	  x)
+ 
+       ;; IFFT by column
+       (scal! (/ (float n 1d0))
+ 	     (dotimes (j (ncols x) x)
+ 	       (declare (type fixnum j))
+ 	       (with-vector-data-addresses ((addr-x (store x))
+ 					    (addr-wsave wsave))
+ 		       (incf-sap :complex-double-float addr-x (* j n))
+ 		       (dfftpack::fortran-zfftb n addr-x addr-wsave)))))))
+ 

Re: [Matlisp-users] FFT/IFFT Seriously Broke for Matrices

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

>>>>> "mak" == mak  <ma...@ll...> writes:

    mak> Okay...so I embarrass myself again...I didn't thoroughly check out the vector
    mak> case, since I was concertrating on the 2D case.  The requires another modification 
    mak> to fft.lisp.  The below DIFF against fft.lisp should be the right one.

    mak> sorry

We should be the ones apologizing!

Patch applied.  I'll commit soon.  Testing with ACL needed, still.

Ray

P.S.  Were you going to send some pseudo inverse code too?

Re: [Matlisp-users] FFT/IFFT Seriously Broke for Matrices

[<ma...@ll...>]

Okay...so I embarrass myself again...I didn't thoroughly check out the vector
case, since I was concertrating on the 2D case.  The requires another modification 
to fft.lisp.  The below DIFF against fft.lisp should be the right one.

sorry

mike



*** src/matlisp/src/fft.lisp	Thu Oct 25 17:48:45 2001
--- src/matlisp-working/src/fft.lisp	Sun Feb  2 12:22:45 2003
***************
*** 217,238 ****
  		       ((col-vector-p x)
  			(make-complex-matrix-dim n 1))
  		       (t (make-complex-matrix-dim n (ncols x))))))
      (if (row-or-col-vector-p x)
  	(progn
  	  (copy! x result)
! 	  (zfftf n (store result) wsave))
  
!       (dotimes (j (ncols x))
! 	(declare (type fixnum j))
! 	 (dotimes (i (nrows x))
! 	   (declare (type fixnum i))
! 	   (setf (matrix-ref result i j) (matrix-ref x i j)))
! 	 (with-vector-data-addresses ((addr-result (store result))
! 				      (addr-wsave wsave))
! 	    (incf-sap :complex-double-float addr-result (* j n))
! 	    (dfftpack::fortran-zfftf n addr-result addr-wsave))))
  
!       result))
  
  
  #+:allegro
--- 217,242 ----
  		       ((col-vector-p x)
  			(make-complex-matrix-dim n 1))
  		       (t (make-complex-matrix-dim n (ncols x))))))
+ 
      (if (row-or-col-vector-p x)
  	(progn
  	  (copy! x result)
! 	  (zfftf n (store result) wsave)
! 	  result)
  
!       ;; Do FFT by column...first we copy the array
!       (progn
! 	;; copy X to a working array, RESULT
! 	(dotimes (j-copy (ncols x))
! 	  (declare (type fixnum j-copy))
! 	  
! 	  ;; Note: we may be asked to truncate if N < (NROWS X)
! 	  (dotimes (i-copy (min n (nrows x))) 
! 	    (declare (type fixnum i-copy))
! 	    (setf (matrix-ref result i-copy j-copy) (matrix-ref x i-copy j-copy))))
  
! 	;; Then do FFT by column
! 	(fft! result)))))
  
  
  #+:allegro
***************
*** 286,313 ****
  	(progn
  	  (copy! x result)
  	  (zfftb n (store result) wsave)
! 	  (let ((scale-factor (/ (float n 1d0))))
! 	    (dotimes (k n)
! 	      (setf (matrix-ref result k)
! 		    (* scale-factor (matrix-ref result k))))))
! 
! 	(let ((scale-factor (/ (float n 1d0))))
! 	  (dotimes (j (ncols x))
! 	    (declare (type fixnum j))
! 	    (dotimes (i (nrows x))
! 	      (declare (type fixnum i))
! 	      (setf (matrix-ref result i j) (matrix-ref x i j)))
! 	    (with-vector-data-addresses ((addr-result (store result))
! 					 (addr-wsave wsave))
! 	      (incf-sap :complex-double-float addr-result (* j n))
! 	      (dfftpack::fortran-zfftb n addr-result addr-wsave))
! 	    ;; Scale the result
! 	    (dotimes (i (nrows x))
! 	      (declare (type fixnum i))
! 	      (setf (matrix-ref result i j) (* scale-factor (matrix-ref x i j)))))))
  
!     result))
  
  
  #+:allegro
  (defmethod ifft ((x standard-matrix) &optional n)
--- 290,310 ----
  	(progn
  	  (copy! x result)
  	  (zfftb n (store result) wsave)
! 	  (scal! (/ (float n 1d0)) result))
  
!       ;; Perform the IFFT by columns...first we copy the array
!       (progn
! 	;; Copy X to the working Array RESULT
! 	(dotimes (j-copy (ncols x))
! 	  (declare (type fixnum j-copy))
! 
! 	  ;; Note: we may be asked to truncate if N < (NROWS X)
! 	  (dotimes (i-copy (min n (nrows x)))
! 	    (declare (type fixnum i-copy))
! 	    (setf (matrix-ref result i-copy j-copy) (matrix-ref x i-copy j-copy))))
  
+ 	;; The we do the IFFT
+ 	(ifft! result)))))
  
  #+:allegro
  (defmethod ifft ((x standard-matrix) &optional n)
***************
*** 348,350 ****
--- 345,406 ----
  
        result))
  
+ ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+ ;; Modification for destructive FFT's.  The purpose is to eliminate
+ ;; cons'ing for certain application requiring numerous calls to these
+ ;; routines.
+ 
+ (defgeneric fft! (x)
+   (:documentation "See FFT but note that the optional N is NOT permitted.
+   Performs in place FFT modifying X.  This will only work for a complex matrix."))
+ 
+ (defgeneric ifft! (x)
+   (:documentation "See IFFT but note that the optional N is NOT permitted.
+   Performs in place IFFT modifying X.  This will only work for a complex matrix."))
+ 
+ #+:cmu  
+ (defmethod fft! ((x complex-matrix))
+   (let* ((n (if (row-or-col-vector-p x)
+ 		(max (nrows x) (ncols x))
+ 	      (nrows x)))
+ 	 (wsave (lookup-wsave-entry n)))
+ 
+     (if (row-or-col-vector-p x)
+ 	(progn
+ 	  (zfftf n (store x) wsave)
+ 	  x)
+ 
+       (dotimes (j (ncols x))
+ 	(declare (type fixnum j))
+ 	(with-vector-data-addresses ((addr-x (store x))
+ 				     (addr-wsave wsave))
+ 	    (incf-sap :complex-double-float addr-x (* j n))
+ 	    (dfftpack::fortran-zfftf n addr-x addr-wsave)))))
+ 
+   ;; return x, the result
+   x)
+ 
+ #+:cmu
+ (defmethod ifft! ((x complex-matrix))
+   (let* ((n (if (row-or-col-vector-p x)
+ 		      (max (nrows x) (ncols x))
+ 		      (nrows x)))
+ 	 (wsave (lookup-wsave-entry n)))
+ 
+     (if (row-or-col-vector-p x)
+ 	(progn
+ 	  (zfftb n (store x) wsave)
+ 	  x)
+ 
+       ;; IFFT by column
+       (dotimes (j (ncols x))
+ 	(declare (type fixnum j))
+ 	(with-vector-data-addresses ((addr-x (store x))
+ 				     (addr-wsave wsave))
+ 	      (incf-sap :complex-double-float addr-x (* j n))
+ 	      (dfftpack::fortran-zfftb n addr-x addr-wsave))))
+ 	    
+     ;; Scale the result
+     (scal! (/ (float n 1d0)) x))
+   x)
+ 

[Matlisp-users] FFT/IFFT Seriously Broke for Matrices

[<ma...@ll...>]

All,

1.  I noted the following today.  FFT & IFFT only works for vectors.
    2D matricies wont compute correctly.  Try (m- (ifft (fft x)) x)
    on x = (ones 16 2).

2.  The problem seems to be a misplaced closing paren on copying the
    input matrix, X, to a working matrix RESULT.

3.  I have made the following change to the #+CMU portions of the code
    (the diff's are attached)

    a.  Add two new functions, FFT! and IFFT!, and exported them.  These
	are useful themselves for performing fast convolution without
	the need for unneeded cons'ing.  I.e. (setf x (ifft (m.* (fft x) y)))
        can be replaced with (ifft! (m.*! y (fft! x)))

    b.  Removed the loops for scaling in IFFT (too much cons'ing via MATRIX-REF)
        Performed scaling with SCAL!...much faster

    c.  FFT and IFFT call FFT! and IFFT! after create copy of matrix.

4.  Can someone look into making similar mods to the allegro portions of the code?

5.  I recommend updating CVS with the changes soon.

6.  Files modified:  FFT.LISP, DFFTPACK.LISP, PACKAGES.LISP

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;		   Changes to PACKAGES.LISP follow
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
*** src/matlisp/packages.lisp	Tue Oct  8 19:30:39 2002
--- src/matlisp-working/packages.lisp	Sat Feb  1 08:02:37 2003
***************
*** 207,213 ****
       "EIG"
       "EYE"
       "FFT"
!      "FFT"
       "FILL-MATRIX"
       "FLOAT-MATRIX"
       "FLOAT-MATRIX-ARRAY-TYPE"
--- 207,213 ----
       "EIG"
       "EYE"
       "FFT"
!      "FFT!"
       "FILL-MATRIX"
       "FLOAT-MATRIX"
       "FLOAT-MATRIX-ARRAY-TYPE"
***************
*** 224,229 ****
--- 224,230 ----
       "GETRS!"
       "HELP"
       "IFFT"
+      "IFFT!"
       "IMAG"
       "JOIN"
       "LOAD-BLAS-&-LAPACK-BINARIES"
***************
*** 285,290 ****
--- 286,292 ----
       "NUMBER-OF-ROWS"
       "ONES"
       "PRINT-ELEMENT"
+      "PINV"
       "QR"
       "QR!"
       "GEQR!"


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;		   Changes to DFFTPACK.LISP follow
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
*** src/matlisp/src/dfftpack.lisp	Tue Jul 11 14:02:03 2000
--- src/matlisp-working/src/dfftpack.lisp	Sat Feb  1 06:07:42 2003
***************
*** 126,132 ****
          destroyed between calls of subroutine zfftf or zfftb
   "        
    (n :integer :input)
!   (c (* :complex-double-float) :output)
    (wsave (* :double-float) :input))
  
  
--- 126,132 ----
          destroyed between calls of subroutine zfftf or zfftb
   "        
    (n :integer :input)
!   (c (* :complex-double-float) :input-output)
    (wsave (* :double-float) :input))
  
  
***************
*** 174,178 ****
          destroyed between calls of subroutine zfftf or zfftb
   "        
    (n :integer :input)
!   (c (* :complex-double-float) :output)
    (wsave (* :double-float) :input))
--- 174,178 ----
          destroyed between calls of subroutine zfftf or zfftb
   "        
    (n :integer :input)
!   (c (* :complex-double-float) :input-output)
    (wsave (* :double-float) :input))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;		   Changes to FFT.LISP follow
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
*** src/matlisp/src/fft.lisp	Thu Oct 25 17:48:45 2001
--- src/matlisp-working/src/fft.lisp	Sat Feb  1 07:39:51 2003
***************
*** 217,238 ****
  		       ((col-vector-p x)
  			(make-complex-matrix-dim n 1))
  		       (t (make-complex-matrix-dim n (ncols x))))))
      (if (row-or-col-vector-p x)
  	(progn
  	  (copy! x result)
  	  (zfftf n (store result) wsave))
  
!       (dotimes (j (ncols x))
! 	(declare (type fixnum j))
! 	 (dotimes (i (nrows x))
! 	   (declare (type fixnum i))
! 	   (setf (matrix-ref result i j) (matrix-ref x i j)))
! 	 (with-vector-data-addresses ((addr-result (store result))
! 				      (addr-wsave wsave))
! 	    (incf-sap :complex-double-float addr-result (* j n))
! 	    (dfftpack::fortran-zfftf n addr-result addr-wsave))))
  
!       result))
  
  
  #+:allegro
--- 217,241 ----
  		       ((col-vector-p x)
  			(make-complex-matrix-dim n 1))
  		       (t (make-complex-matrix-dim n (ncols x))))))
+ 
      (if (row-or-col-vector-p x)
  	(progn
  	  (copy! x result)
  	  (zfftf n (store result) wsave))
  
!       ;; Do FFT by column...first we copy the array
!       (progn
! 	;; copy X to a working array, RESULT
! 	(dotimes (j-copy (ncols x))
! 	  (declare (type fixnum j-copy))
! 	  
! 	  ;; Note: we may be asked to truncate if N < (NROWS X)
! 	  (dotimes (i-copy (min n (nrows x))) 
! 	    (declare (type fixnum i-copy))
! 	    (setf (matrix-ref result i-copy j-copy) (matrix-ref x i-copy j-copy))))
  
! 	;; Then do FFT by column
! 	(fft! result)))))
  
  
  #+:allegro
***************
*** 286,313 ****
  	(progn
  	  (copy! x result)
  	  (zfftb n (store result) wsave)
! 	  (let ((scale-factor (/ (float n 1d0))))
! 	    (dotimes (k n)
! 	      (setf (matrix-ref result k)
! 		    (* scale-factor (matrix-ref result k))))))
! 
! 	(let ((scale-factor (/ (float n 1d0))))
! 	  (dotimes (j (ncols x))
! 	    (declare (type fixnum j))
! 	    (dotimes (i (nrows x))
! 	      (declare (type fixnum i))
! 	      (setf (matrix-ref result i j) (matrix-ref x i j)))
! 	    (with-vector-data-addresses ((addr-result (store result))
! 					 (addr-wsave wsave))
! 	      (incf-sap :complex-double-float addr-result (* j n))
! 	      (dfftpack::fortran-zfftb n addr-result addr-wsave))
! 	    ;; Scale the result
! 	    (dotimes (i (nrows x))
! 	      (declare (type fixnum i))
! 	      (setf (matrix-ref result i j) (* scale-factor (matrix-ref x i j)))))))
  
!     result))
  
  
  #+:allegro
  (defmethod ifft ((x standard-matrix) &optional n)
--- 289,309 ----
  	(progn
  	  (copy! x result)
  	  (zfftb n (store result) wsave)
! 	  (scal! (/ (float n 1d0)) result))
  
!       ;; Perform the IFFT by columns...first we copy the array
!       (progn
! 	;; Copy X to the working Array RESULT
! 	(dotimes (j-copy (ncols x))
! 	  (declare (type fixnum j-copy))
! 
! 	  ;; Note: we may be asked to truncate if N < (NROWS X)
! 	  (dotimes (i-copy (min n (nrows x)))
! 	    (declare (type fixnum i-copy))
! 	    (setf (matrix-ref result i-copy j-copy) (matrix-ref x i-copy j-copy))))
  
+ 	;; The we do the IFFT
+ 	(ifft! result)))))
  
  #+:allegro
  (defmethod ifft ((x standard-matrix) &optional n)
***************
*** 348,350 ****
--- 344,403 ----
  
        result))
  
+ ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+ ;; Modification for destructive FFT's.  The purpose is to eliminate
+ ;; cons'ing for certain application requiring numerous calls to these
+ ;; routines.
+ 
+ (defgeneric fft! (x)
+   (:documentation "See FFT but note that the optional N is NOT permitted.
+   Performs in place FFT modifying X.  This will only work for a complex matrix."))
+ 
+ (defgeneric ifft! (x)
+   (:documentation "See IFFT but note that the optional N is NOT permitted.
+   Performs in place IFFT modifying X.  This will only work for a complex matrix."))
+ 
+ #+:cmu  
+ (defmethod fft! ((x complex-matrix))
+   (let* ((n (if (row-or-col-vector-p x)
+ 		(max (nrows x) (ncols x))
+ 	      (nrows x)))
+ 	 (wsave (lookup-wsave-entry n)))
+ 
+     (if (row-or-col-vector-p x)
+ 	(zfftf n (store x) wsave)
+ 
+       (dotimes (j (ncols x))
+ 	(declare (type fixnum j))
+ 	(with-vector-data-addresses ((addr-x (store x))
+ 				     (addr-wsave wsave))
+ 	    (incf-sap :complex-double-float addr-x (* j n))
+ 	    (dfftpack::fortran-zfftf n addr-x addr-wsave)))))
+ 
+   ;; return x, the result
+   x)
+ 
+ #+:cmu
+ (defmethod ifft! ((x complex-matrix))
+   (let* ((n (if (row-or-col-vector-p x)
+ 		      (max (nrows x) (ncols x))
+ 		      (nrows x)))
+ 	 (wsave (lookup-wsave-entry n)))
+ 
+     (if (row-or-col-vector-p x)
+ 	(progn
+ 	  (zfftb n (store x) wsave)
+ 	  (scal! (/ (float n 1d0)) x))
+ 
+       ;; IFFT by column
+       (dotimes (j (ncols x))
+ 	(declare (type fixnum j))
+ 	(with-vector-data-addresses ((addr-x (store x))
+ 				     (addr-wsave wsave))
+ 	      (incf-sap :complex-double-float addr-x (* j n))
+ 	      (dfftpack::fortran-zfftb n addr-x addr-wsave))))
+ 	    
+     ;; Scale the result
+     (scal! (/ (float n 1d0)) x))
+   x)
+ 

Re: [Matlisp-users] PINV, initial release

[Jefferson P. <jp...@cs...>]

Raymond Toy wrote:

> 
> o Should the args to pinv be keyword args instead of optional?  If so,
>   please make the mthd arg have the keyword :method instead.  Easier
>   to figure out what it is. :-)

I thought that the MATLISP convention was to use optional arguments. 
This matches SVD, JOIN, etc.  Personally, I almost always prefer keyword 
arguments over optional.  Especially when there's more than one optional 
arg.

I have no real preference on the name, between PNV and pseudoinverse. 
PINV is nice because it matches matlab/octave, which makes it easy to 
find, if you're looking for it.  But pseudoinverse will be equally easy, 
I think.

J.

Re: [Matlisp-users] PINV, initial release

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

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

    Tunc> Hi Mike;  I can enter the changes.  But if you're going
    Tunc> to update PINV in few weeks I'd rather do it then, to avoid
    Tunc> doing it twice.  --Thanks, Tunc

Sounds good to me.

A couple of comments:

o Can this be called pseudo-inverse instead of pinv?  My general rule
  has been if LAPACK has a certain useful routine, then use that name,
  otherwise pick a more useful, less abbreviated name. :-)

o Should the args to pinv be keyword args instead of optional?  If so,
  please make the mthd arg have the keyword :method instead.  Easier
  to figure out what it is. :-)

Everything else looks good.

Ray

Re: [Matlisp-users] PINV, initial release

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

Hi Mike;  I can enter the changes.  But if you're going
to update PINV in few weeks I'd rather do it then, to avoid
doing it twice.  --Thanks, Tunc

ma...@ll... wrote:
> 
> All,
> 
> 1.  Below are proposed changes to add pseudo-inverse, PINV, to matlisp.
> 
> 2.  I've tested this on some random real and complex matrix and compared
>     the results to OCTAVE including rank deficient matricies (they match).
> 
> 3.  My intentions are to replace PINV in the next few weeks with direct
>     calls to DGELSD.F/ZGELSD.F.  I intend _no_ change to the calling or
>     return arguments.
> 
> 4.  I'd specifically like feedback on the inclusion of the
>     LEAST-SQUARES option; MTHD.  I almost removed this block of code
>     from PINV several times as being too trivial and adding, otherwise,
>     unnecessary code.  Its in the below code only because I posted a
>     similar routine last week.
> 
> 5.  Ray, Tunc: If there are no objections, or agreed to changes I need
>     to make, by the end of the week of 12/2/02, could you enter the
>     below changes in CVS?
> 
> tnx
> 
> mike
> 
> ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
> ;;;                     CHANGES TO SYSTEM.DCL
> ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
> diff -c ../matlisp system.dcl
> 
> *** ../matlisp/system.dcl       Tue Oct  8 19:30:39 2002
> --- system.dcl  Fri Nov 29 11:39:53 2002
> ***************
> *** 218,224 ****
>                      "mdivide"
>                      "msqrt"
>                      #-:mswindows "fft"
> !                    "geqr"))
>          (:module "special-functions"
>                 :source-pathname "matlisp:src"
>                 :binary-pathname ""
> --- 218,225 ----
>                      "mdivide"
>                      "msqrt"
>                      #-:mswindows "fft"
> !                    "geqr"
> !                    "pinv"))
>          (:module "special-functions"
>                 :source-pathname "matlisp:src"
>                 :binary-pathname ""
> 
> ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
> ;;;                    CHANGES TO PACKAGES.LISP
> ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
> diff -c ../matlisp packages.lisp
> 
> *** ../matlisp/packages.lisp    Tue Oct  8 19:30:39 2002
> --- packages.lisp       Fri Nov 29 11:36:07 2002
> ***************
> *** 285,290 ****
> --- 285,291 ----
>        "NUMBER-OF-ROWS"
>        "ONES"
>        "PRINT-ELEMENT"
> +      "PINV"
>        "QR"
>        "QR!"
>        "GEQR!"
> 
> ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
> ;;;               NEW FILE FOR matlisp/src/pinv.lisp
> ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
> (in-package "MATLISP")
> 
> ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
> (defun pinv (a &optional (b nil) (tol nil) (mthd :svd))
>   "
>   SYNTAX
>   ======
>   (PINV A [B TOL MTHD])
> 
>   INPUT
>   -----
>   A         A Matlisp matrix, M x N
>   B         Optional.  Default NIL.  If provided, return solution to
>             A X = B.
>   TOL       Optional.  Default NIL.  Set limit on range of singular values
>             used in computation of pseudo-inverse.
>   MTHD      Optional.  Default :SVD.  Indicates the type of computation
>             to use for the pseudo-inverse.
> 
>               :LS  Force least squares interpretation
>               :SVD Force SVD based computation (default)
> 
>   OUTPUT  (VALUES PINV RANK SIGMA)
>   ------
>   PINV      Pseudo-inverse or solution to A X = B, a Matlisp matrix.
>   RANK      The effective rank of A, r = rank(A).  NIL if MTHD is :LS
>   SIGMA     The list of lenght min(M,N) of the signular values of A.
>             NIL if MTHD is :LS
> 
>   DESCRIPTION
>   ===========
>   Given the equation A X = B, solve for the value of A^+ which will
>   solve X = A^+ B.  A is assumed M x N.  X is N x L and B M x L.  The
>   solution is computed via a signular value decomposition, SVD.  Let
> 
>   A = U S V^H
> 
>   define S^+ = diag(1/s1, 1/s2, ... 1/sr, 0 ... 0) where r is rank(A)
>   Then the p-inv, A^+ = V S^+ U^H.
> 
>   The computuation of r = rank(A) is controlled by TOL.  If TOL is NIL,
>   or TOL < 0, or TOL > 1, TOL is replaced with DOUBLE-FLOAT-EPSILON.
>   All signular values <= TOL*S(1) are considered to be equal
>   to zero.  (Note: singular values have a descending order
>   S(1) >= S(2) ...S(r) > 0.)
> 
>   If A is M x M, and full rank, (M:M/ A B) can be used.
> 
>   Use MTHD :SVD when rank(A) = r < N.  For r <= M < N, this is the
>   underconstrained case for which A^+ is the minimum norm solution.
> 
>   The MTHD :LS can be used when M > N and rank(A) = N.  This is the
>   canonical least squares problem.  In this case A^+ = inv(A^H A) A^H.
>   NOTE:  NO RANK CHECKS OR SUCH ARE MADE HERE.  IF A^H A IS RANK
>   DEFICIENT THIS METHOD MAY FAIL.  USER BEWARE."
> 
>   ;; reset the TOL value as required.  This odd screening
>   ;; is for future compatibility to the ?GELSD.F routines
>   (if (or (null tol)
>           (>= tol 1)
>           (< tol 0))
>       (setf tol double-float-epsilon))
> 
>   (cond
>     ;; Force the least squares interpretation.
>    ((eq mthd :ls)
>     (let* ((ah (ctranspose a))
>            (pinv (m/ (m* ah a) ah)))
>       (values (if b (m* pinv b) pinv) nil nil)))
> 
>    ;; This is the major routine here, using the SVD method
>    ((eq mthd :svd)
>     (multiple-value-bind(up sp vptranspose info)
>         (svd a :s)                      ; use economy SVD
> 
>       (if (null info) (error "Computation of SVD failed."))
> 
>       ;; This form sets all signular values above the rank
>       ;; of A to zero
>       (let ((pinv nil)
>             (sp-values (m::store (diag sp)))
>             (rank
>              ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
>              ;; This form computes the reciprocal of the singular values
>              ;; up to and including the rank, r, and returns the RANK
>              (do ((limit (* tol (matrix-ref sp 0 0)))
>                   (max-rank (min (number-of-rows sp) (number-of-cols sp)))
>                   (n 0 (incf n)))
>                  ;; terminate when all sigular values are used,
>                  ;; or fall below the limit
>                  ((or (= n max-rank)
>                       (<= (matrix-ref sp n n) limit)) n)
> 
>                (setf (matrix-ref sp n n) (/ (matrix-ref sp n n))))))
> 
>         ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
>         ;; set the values of the reciprocal singular values to zero
>         ;; above the RANK of A
>         (do ((n rank (incf n))
>              (max-rank (min (number-of-rows sp) (number-of-cols sp))))
>             ((= n max-rank))
> 
>           (setf (matrix-ref sp n n) 0d0))
> 
>         ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
>         ;; compute the PINV or PINV * B
>         (setf pinv (m* (ctranspose vptranspose) (m* sp (ctranspose up))))
>         (values
>          (if b (m* pinv b) pinv)
>          rank
>          sp-values))))
> 
>    (t (error "Invalid Method (~A) passed to PINV." mthd))))
> 
> -------------------------------------------------------
> This SF.net email is sponsored by: Get the new Palm Tungsten T
> handheld. Power & Color in a compact size!
> http://ads.sourceforge.net/cgi-bin/redirect.pl?palm0002en
> _______________________________________________
> Matlisp-users mailing list
> Mat...@li...
> https://lists.sourceforge.net/lists/listinfo/matlisp-users

[Matlisp-users] PINV, initial release

[<ma...@ll...>]

All,

1.  Below are proposed changes to add pseudo-inverse, PINV, to matlisp.

2.  I've tested this on some random real and complex matrix and compared
    the results to OCTAVE including rank deficient matricies (they match).

3.  My intentions are to replace PINV in the next few weeks with direct
    calls to DGELSD.F/ZGELSD.F.  I intend _no_ change to the calling or
    return arguments.

4.  I'd specifically like feedback on the inclusion of the
    LEAST-SQUARES option; MTHD.  I almost removed this block of code
    from PINV several times as being too trivial and adding, otherwise,
    unnecessary code.  Its in the below code only because I posted a
    similar routine last week.

5.  Ray, Tunc: If there are no objections, or agreed to changes I need
    to make, by the end of the week of 12/2/02, could you enter the
    below changes in CVS?

tnx

mike


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;			CHANGES TO SYSTEM.DCL
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
diff -c ../matlisp system.dcl 

*** ../matlisp/system.dcl	Tue Oct  8 19:30:39 2002
--- system.dcl	Fri Nov 29 11:39:53 2002
***************
*** 218,224 ****
  		     "mdivide"
  		     "msqrt"
  		     #-:mswindows "fft"
! 		     "geqr"))
         (:module "special-functions"
  		:source-pathname "matlisp:src"
  		:binary-pathname ""
--- 218,225 ----
  		     "mdivide"
  		     "msqrt"
  		     #-:mswindows "fft"
! 		     "geqr"
! 		     "pinv"))
         (:module "special-functions"
  		:source-pathname "matlisp:src"
  		:binary-pathname ""

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;		       CHANGES TO PACKAGES.LISP
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
diff -c ../matlisp packages.lisp

*** ../matlisp/packages.lisp	Tue Oct  8 19:30:39 2002
--- packages.lisp	Fri Nov 29 11:36:07 2002
***************
*** 285,290 ****
--- 285,291 ----
       "NUMBER-OF-ROWS"
       "ONES"
       "PRINT-ELEMENT"
+      "PINV"
       "QR"
       "QR!"
       "GEQR!"

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;		  NEW FILE FOR matlisp/src/pinv.lisp
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(in-package "MATLISP")

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun pinv (a &optional (b nil) (tol nil) (mthd :svd))
  "
  SYNTAX
  ======
  (PINV A [B TOL MTHD])

  INPUT
  -----
  A         A Matlisp matrix, M x N
  B         Optional.  Default NIL.  If provided, return solution to
            A X = B.
  TOL       Optional.  Default NIL.  Set limit on range of singular values
            used in computation of pseudo-inverse.
  MTHD      Optional.  Default :SVD.  Indicates the type of computation
            to use for the pseudo-inverse.  

              :LS  Force least squares interpretation
              :SVD Force SVD based computation (default)

  OUTPUT  (VALUES PINV RANK SIGMA)
  ------
  PINV      Pseudo-inverse or solution to A X = B, a Matlisp matrix.
  RANK      The effective rank of A, r = rank(A).  NIL if MTHD is :LS
  SIGMA     The list of lenght min(M,N) of the signular values of A.
            NIL if MTHD is :LS

  DESCRIPTION
  ===========
  Given the equation A X = B, solve for the value of A^+ which will
  solve X = A^+ B.  A is assumed M x N.  X is N x L and B M x L.  The
  solution is computed via a signular value decomposition, SVD.  Let

  A = U S V^H

  define S^+ = diag(1/s1, 1/s2, ... 1/sr, 0 ... 0) where r is rank(A)
  Then the p-inv, A^+ = V S^+ U^H.

  The computuation of r = rank(A) is controlled by TOL.  If TOL is NIL,
  or TOL < 0, or TOL > 1, TOL is replaced with DOUBLE-FLOAT-EPSILON.
  All signular values <= TOL*S(1) are considered to be equal
  to zero.  (Note: singular values have a descending order
  S(1) >= S(2) ...S(r) > 0.)

  If A is M x M, and full rank, (M:M/ A B) can be used.

  Use MTHD :SVD when rank(A) = r < N.  For r <= M < N, this is the
  underconstrained case for which A^+ is the minimum norm solution.

  The MTHD :LS can be used when M > N and rank(A) = N.  This is the
  canonical least squares problem.  In this case A^+ = inv(A^H A) A^H.
  NOTE:  NO RANK CHECKS OR SUCH ARE MADE HERE.  IF A^H A IS RANK
  DEFICIENT THIS METHOD MAY FAIL.  USER BEWARE."

  ;; reset the TOL value as required.  This odd screening
  ;; is for future compatibility to the ?GELSD.F routines
  (if (or (null tol)
	  (>= tol 1)
	  (< tol 0))
      (setf tol double-float-epsilon))

  (cond
    ;; Force the least squares interpretation.
   ((eq mthd :ls)
    (let* ((ah (ctranspose a))
	   (pinv (m/ (m* ah a) ah)))
      (values (if b (m* pinv b) pinv) nil nil)))

   ;; This is the major routine here, using the SVD method
   ((eq mthd :svd)
    (multiple-value-bind(up sp vptranspose info)
	(svd a :s)			; use economy SVD
      
      (if (null info) (error "Computation of SVD failed."))

      ;; This form sets all signular values above the rank
      ;; of A to zero
      (let ((pinv nil)
	    (sp-values (m::store (diag sp)))
	    (rank
	     ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	     ;; This form computes the reciprocal of the singular values 
	     ;; up to and including the rank, r, and returns the RANK
	     (do ((limit (* tol (matrix-ref sp 0 0)))
		  (max-rank (min (number-of-rows sp) (number-of-cols sp)))
		  (n 0 (incf n)))
		 ;; terminate when all sigular values are used,
		 ;; or fall below the limit
		 ((or (= n max-rank)
		      (<= (matrix-ref sp n n) limit)) n)
	       
	       (setf (matrix-ref sp n n) (/ (matrix-ref sp n n))))))

	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	;; set the values of the reciprocal singular values to zero
	;; above the RANK of A
	(do ((n rank (incf n))
	     (max-rank (min (number-of-rows sp) (number-of-cols sp))))
	    ((= n max-rank))

	  (setf (matrix-ref sp n n) 0d0))

	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
	;; compute the PINV or PINV * B 
	(setf pinv (m* (ctranspose vptranspose) (m* sp (ctranspose up))))
	(values
	 (if b (m* pinv b) pinv)
	 rank
	 sp-values))))

   (t (error "Invalid Method (~A) passed to PINV." mthd))))

Re: [Matlisp-users] pseudoinverse

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

Tunc,

> Hi Mike.  This sounds good.  Will the function have to
> work for A,B complex (or A real, B complex and so forth).

Yes...all of those cases will be covered.

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

Re: [Matlisp-users] pseudoinverse

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

Hi Mike.  This sounds good.  Will the function have to
work for A,B complex (or A real, B complex and so forth).
Thanks, Tunc


mak%koe...@ll... wrote:
> 
> Tunc,
> 
> 1.  Over the weekend I spent some more time on the PINV question.  In
> response to your question about DGELSD, the more I thought about it
> the more I liked it.  Since the routine accepts multiple RHS, I
> suppose it can be provided an identity matrix.  The result should then
> be the PINV matrix.  If I write something like this
> 
>   (DEFUN PINV (A &OPTIONAL (B NIL) (...))
>      ...
>      (LET ((B (IF (NULL B) IDENTITY-FOR-A B)))
>        ...
>        (CALLING-GELSD ...) ...))
> 
> then the user can say (PINV A) to get the psuedo-inverse,
> or they can say (PINV A B) to solve the multiple RHS problem.
> 
> 2.  The use of ?GELSD should aviod (my) concerns over differences in
> D1MACH values in a given CL implimentation and the associated FORTRAN.
> (BTW I found that the machine values are controlled by DLAMCH.F for
> LAPACK, not D1MACH.F as I thought.)
> 
> 3.  Due to current time constraints, I intend to:
> 
>     a. change the calling convention of PINV I posted last week to
>        provide this capability.  The computations remain essentially
>        unchanged.
> 
>     b. over the next few weeks, develop the ?GELSD.F wrappers and
>        replace the internals of PINV to use these.
> 
>     any thoughts or redirections concerning this?
> 
> mike
> 
> -------------------------------------------------------
> 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] pseudoinverse

[mak%<koe...@ll...>]

Tunc,

1.  Over the weekend I spent some more time on the PINV question.  In
response to your question about DGELSD, the more I thought about it
the more I liked it.  Since the routine accepts multiple RHS, I
suppose it can be provided an identity matrix.  The result should then
be the PINV matrix.  If I write something like this

  (DEFUN PINV (A &OPTIONAL (B NIL) (...))
     ...
     (LET ((B (IF (NULL B) IDENTITY-FOR-A B)))
       ...
       (CALLING-GELSD ...) ...))

then the user can say (PINV A) to get the psuedo-inverse,
or they can say (PINV A B) to solve the multiple RHS problem.

2.  The use of ?GELSD should aviod (my) concerns over differences in
D1MACH values in a given CL implimentation and the associated FORTRAN.
(BTW I found that the machine values are controlled by DLAMCH.F for
LAPACK, not D1MACH.F as I thought.)

3.  Due to current time constraints, I intend to:

    a. change the calling convention of PINV I posted last week to
       provide this capability.  The computations remain essentially
       unchanged.

    b. over the next few weeks, develop the ?GELSD.F wrappers and 
       replace the internals of PINV to use these.


    any thoughts or redirections concerning this?

mike

Re: [Matlisp-users] pseudoinverse

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

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

    Michael> I haven't studied this enough yet, but it appears that D1MACH.LISP is
    Michael> used via F2CL to provide MACH information to the FORTRAN routines in,
    Michael> e.g., TOMS.  D1MACH draws its information from the CL implimentation
    Michael> limits.  I don't know that this is the same as that produced by
    Michael> underlying FORTRAN compiler's treatement of D1MACH.F.  Unless it
    Michael> is the same, the numerically accuracy of the routines that rely
    Michael> on it could be impacted.  This is just a concern on my part, not 

Yes, this is true in general, but if LAPACK needs stuff from d1mach,
it's going to get it from the Fortran d1mach, not the Lisp d1mach.

There is a problem, however, because if there is such a thing for
LAPACK, it probably hasn't been tweaked at all so whatever those the
defaults are are getting used.

Ray

Re: [Matlisp-users] pseudoinverse

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

Tunc,


> Can you say something about the difference between
> your approach (using SVD) and the LAPACK routines called
> DGELS and DGELSD. Numerical stability?

DGELS will handle only full rank A, DGELSD will solve multiple
right hand sides of Ax = [b1 | b2 | ... ] (each bi a column vector).
An entity called PINV is not returned.  If PINV is not required,
then DGELSD might be a better choice.

> 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). 

I haven't studied this enough yet, but it appears that D1MACH.LISP is
used via F2CL to provide MACH information to the FORTRAN routines in,
e.g., TOMS.  D1MACH draws its information from the CL implimentation
limits.  I don't know that this is the same as that produced by
underlying FORTRAN compiler's treatement of D1MACH.F.  Unless it
is the same, the numerically accuracy of the routines that rely
on it could be impacted.  This is just a concern on my part, not 
something that I have even looked at yet.  Perhaps you or Ray can
comment.

mike

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] 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

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

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

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

If you want to know how to extend matlisp, please ask and ye shall
receive.

What do you want to discuss about add-ons?  How to do them?  Whether
they should be added?  Useful additions?

Ray

Re: [Matlisp-users] pseudoinverse

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

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

> Tunc,
> 
> > When M<N, should the matrix have rank=M or does it
> > not matter?

Tunc,

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
for users to learn how matlisp is extended.

Nicolas.

Re: [Matlisp-users] pseudoinverse

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

Tunc,

> When M<N, should the matrix have rank=M or does it
> not matter?

1.  Good question.  The answer is, PINV as written will fail or
    produce a numerically unstable result due to taking the reciprocal
    of a very small singular value.

2.  The solution is straight forward, let all singular values < MACH-TOL
    equal zero.  ==> rank(A) = r <= min(M,N).  And use just those 
    r elements to form PINV.

3.  I have an all day meeting today, but will make these changes in the
    next couple of days.  ....BTW what is the name of the "*MACH" in 
    the matlisp core....

more later


mike