[lisplab-cvs] r116 - src/matlisp src/matrix
Jørn Inge Vestgården
jivestgarden at common-lisp.net
Sun Nov 15 18:11:20 UTC 2009
Author: jivestgarden
Date: Sun Nov 15 13:11:19 2009
New Revision: 116
Log:
added tridag lin-solve
Added:
src/matlisp/tridiag.lisp
Modified:
lisplab.asd
package.lisp
src/matlisp/lapack.lisp
src/matrix/level1-dgt.lisp
Modified: lisplab.asd
==============================================================================
--- lisplab.asd (original)
+++ lisplab.asd Sun Nov 15 13:11:19 2009
@@ -178,7 +178,8 @@
(:file "lapack")
(:file "mul")
(:file "inv")
- (:file "geev")))))
+ (:file "geev")
+ (:file "tridiag")))))
(defsystem :lisplab-fftw
:depends-on (:lisplab-base)
Modified: package.lisp
==============================================================================
--- package.lisp (original)
+++ package.lisp Sun Nov 15 13:11:19 2009
@@ -133,6 +133,8 @@
"MATRIX-GE"
"MATRIX-DGE"
"MATRIX-ZGE"
+ "MATRIX-DGT"
+ "MATRIX-DDI"
"FUNCTION-MATRIX"
"MATRIX-SPARSE"
"*LISPLAB-PRINT-SIZE*"
Modified: src/matlisp/lapack.lisp
==============================================================================
--- src/matlisp/lapack.lisp (original)
+++ src/matlisp/lapack.lisp Sun Nov 15 13:11:19 2009
@@ -1628,4 +1628,243 @@
(IPIV (* :integer) :input)
(WORK (* :complex-double-float) :workspace-output)
(LWORK :integer :input)
- (INFO :integer :output))
\ No newline at end of file
+ (INFO :integer :output))
+
+(def-fortran-routine dgttrf :void
+ "-- LAPACK routine (version 3.2) --
+ -- LAPACK is a software package provided by Univ. of Tennessee, --
+ -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+ November 2006
+
+ .. Scalar Arguments ..
+ INTEGER INFO, N
+ ..
+ .. Array Arguments ..
+ INTEGER IPIV( * )
+ DOUBLE PRECISION D( * ), DL( * ), DU( * ), DU2( * )
+ ..
+
+ Purpose
+ =======
+
+ DGTTRF computes an LU factorization of a real tridiagonal matrix A
+ using elimination with partial pivoting and row interchanges.
+
+ The factorization has the form
+ A = L * U
+ where L is a product of permutation and unit lower bidiagonal
+ matrices and U is upper triangular with nonzeros in only the main
+ diagonal and first two superdiagonals.
+
+ Arguments
+ =========
+
+ N (input) INTEGER
+ The order of the matrix A.
+
+ DL (input/output) DOUBLE PRECISION array, dimension (N-1)
+ On entry, DL must contain the (n-1) sub-diagonal elements of
+ A.
+
+ On exit, DL is overwritten by the (n-1) multipliers that
+ define the matrix L from the LU factorization of A.
+
+ D (input/output) DOUBLE PRECISION array, dimension (N)
+ On entry, D must contain the diagonal elements of A.
+
+ On exit, D is overwritten by the n diagonal elements of the
+ upper triangular matrix U from the LU factorization of A.
+
+ DU (input/output) DOUBLE PRECISION array, dimension (N-1)
+ On entry, DU must contain the (n-1) super-diagonal elements
+ of A.
+
+ On exit, DU is overwritten by the (n-1) elements of the first
+ super-diagonal of U.
+
+ DU2 (output) DOUBLE PRECISION array, dimension (N-2)
+ On exit, DU2 is overwritten by the (n-2) elements of the
+ second super-diagonal of U.
+
+ IPIV (output) INTEGER array, dimension (N)
+ The pivot indices; for 1 <= i <= n, row i of the matrix was
+ interchanged with row IPIV(i). IPIV(i) will always be either
+ i or i+1; IPIV(i) = i indicates a row interchange was not
+ required.
+
+ INFO (output) INTEGER
+ = 0: successful exit
+ < 0: if INFO = -k, the k-th argument had an illegal value
+ > 0: if INFO = k, U(k,k) is exactly zero. The factorization
+ has been completed, but the factor U is exactly
+ singular, and division by zero will occur if it is used
+ to solve a system of equations."
+ (N :integer :input)
+ (DL (* :double-float) :input-output)
+ (D (* :double-float) :input-output)
+ (DU (* :double-float) :input-output)
+ (DU2 (* :double-float) :input-output)
+ (IPIV (* :integer) :input)
+ (INFO :integer :output))
+
+(def-fortran-routine dgttrs :void
+ "-- LAPACK routine (version 3.2) --
+ -- LAPACK is a software package provided by Univ. of Tennessee, --
+ -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+ November 2006
+
+ .. Scalar Arguments ..
+ CHARACTER TRANS
+ INTEGER INFO, LDB, N, NRHS
+ ..
+ .. Array Arguments ..
+ INTEGER IPIV( * )
+ DOUBLE PRECISION B( LDB, * ), D( * ), DL( * ), DU( * ), DU2( * )
+ ..
+
+ Purpose
+ =======
+
+ DGTTRS solves one of the systems of equations
+ A*X = B or A'*X = B,
+ with a tridiagonal matrix A using the LU factorization computed
+ by DGTTRF.
+
+ 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.
+
+ NRHS (input) INTEGER
+ The number of right hand sides, i.e., the number of columns
+ of the matrix B. NRHS >= 0.
+
+ DL (input) DOUBLE PRECISION array, dimension (N-1)
+ The (n-1) multipliers that define the matrix L from the
+ LU factorization of A.
+
+ D (input) DOUBLE PRECISION array, dimension (N)
+ The n diagonal elements of the upper triangular matrix U from
+ the LU factorization of A.
+
+ DU (input) DOUBLE PRECISION array, dimension (N-1)
+ The (n-1) elements of the first super-diagonal of U.
+
+ DU2 (input) DOUBLE PRECISION array, dimension (N-2)
+ The (n-2) elements of the second super-diagonal of U.
+
+ IPIV (input) INTEGER array, dimension (N)
+ The pivot indices; for 1 <= i <= n, row i of the matrix was
+ interchanged with row IPIV(i). IPIV(i) will always be either
+ i or i+1; IPIV(i) = i indicates a row interchange was not
+ required.
+
+ B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
+ On entry, the matrix of right hand side vectors B.
+ On exit, B is overwritten by the solution vectors 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)
+ (NRHS :integer :input)
+ (DL (* :double-float) :input)
+ (D (* :double-float) :input)
+ (DU (* :double-float) :input)
+ (DU2 (* :double-float) :input)
+ (IPIV (* :integer) :input)
+ (B (* :double-float) :input-output)
+ (LDB :integer :input)
+ (INFO :integer :output))
+
+(def-fortran-routine dgtsv :void
+ "-- LAPACK routine (version 3.2) --
+ -- LAPACK is a software package provided by Univ. of Tennessee, --
+ -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+ November 2006
+
+ .. Scalar Arguments ..
+ INTEGER INFO, LDB, N, NRHS
+ ..
+ .. Array Arguments ..
+ DOUBLE PRECISION B( LDB, * ), D( * ), DL( * ), DU( * )
+ ..
+
+ Purpose
+ =======
+
+ DGTSV solves the equation
+
+ A*X = B,
+
+ where A is an n by n tridiagonal matrix, by Gaussian elimination with
+ partial pivoting.
+
+ Note that the equation A'*X = B may be solved by interchanging the
+ order of the arguments DU and DL.
+
+ Arguments
+ =========
+
+ 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.
+
+ DL (input/output) DOUBLE PRECISION array, dimension (N-1)
+ On entry, DL must contain the (n-1) sub-diagonal elements of
+ A.
+
+ On exit, DL is overwritten by the (n-2) elements of the
+ second super-diagonal of the upper triangular matrix U from
+ the LU factorization of A, in DL(1), ..., DL(n-2).
+
+ D (input/output) DOUBLE PRECISION array, dimension (N)
+ On entry, D must contain the diagonal elements of A.
+
+ On exit, D is overwritten by the n diagonal elements of U.
+
+ DU (input/output) DOUBLE PRECISION array, dimension (N-1)
+ On entry, DU must contain the (n-1) super-diagonal elements
+ of A.
+
+ On exit, DU is overwritten by the (n-1) elements of the first
+ super-diagonal of U.
+
+ B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
+ On entry, the N by NRHS matrix of right hand side matrix B.
+ On exit, if INFO = 0, the N by NRHS 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
+ > 0: if INFO = i, U(i,i) is exactly zero, and the solution
+ has not been computed. The factorization has not been
+ completed unless i = N."
+ (N :integer :input)
+ (NRHS :integer :input)
+ (DL (* :double-float) :input-output)
+ (D (* :double-float) :input-output)
+ (DU (* :double-float) :input-output)
+ (B (* :double-float) :input-output)
+ (LDB :integer :input)
+ (INFO :integer :output))
+
+
+
Added: src/matlisp/tridiag.lisp
==============================================================================
--- (empty file)
+++ src/matlisp/tridiag.lisp Sun Nov 15 13:11:19 2009
@@ -0,0 +1,38 @@
+;;; Lisplab, matliap/tridiag.lisp
+;;; Lapack-based, tridiagonal routines
+
+;;; Copyright (C) 2009 Joern Inge Vestgaarden
+;;;
+;;; This program is free software; you can redistribute it and/or modify
+;;; it under the terms of the GNU General Public License as published by
+;;; the Free Software Foundation; either version 2 of the License, or
+;;; (at your option) any later version.
+;;;
+;;; This program is distributed in the hope that it will be useful,
+;;; but WITHOUT ANY WARRANTY; without even the implied warranty of
+;;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+;;; GNU General Public License for more details.
+;;;
+;;; You should have received a copy of the GNU General Public License along
+;;; with this program; if not, write to the Free Software Foundation, Inc.,
+;;; 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+(in-package :lisplab)
+
+(defmethod lin-solve ((A matrix-lisp-dgt) (b matrix-blas-dge))
+ (lin-solve! (copy A) (copy b)))
+
+(defmethod lin-solve! ((A matrix-lisp-dgt) (b matrix-blas-dge))
+ ;; TODO catch error from input
+ (if cl-user::*lisplab-liblapack-path*
+ (let* ((N (cols A)))
+ (f77::dgtsv N
+ 1
+ (slot-value a 'subdiagonal-store )
+ (slot-value a 'diagonal-store)
+ (slot-value a 'superdiagonal-store)
+ (matrix-store b)
+ N
+ 0)
+ b)
+ (call-next-method)))
\ No newline at end of file
Modified: src/matrix/level1-dgt.lisp
==============================================================================
--- src/matrix/level1-dgt.lisp (original)
+++ src/matrix/level1-dgt.lisp Sun Nov 15 13:11:19 2009
@@ -95,8 +95,7 @@
((= (1+ row) col)
(setf (aref (slot-value matrix 'superdiagonal-store) row)
val2))
- (t
- (warn "Array out of bonds for tridiagonal matrix. Ignored.")))))
+ #+nil (t (warn "Array out of bonds for tridiagonal matrix. Ignored.")))))
(defmethod vref ((matrix matrix-base-dgt) idx)
(let ((len (slot-value matrix 'rowcols)))
@@ -106,7 +105,7 @@
(aref (slot-value matrix 'superdiagonal-store) (- idx len)))
((< idx (slot-value matrix 'size))
(aref (slot-value matrix 'subdiagonal-store) (- idx (- (* 2 len) 1))))
- (t
+ #+nil (t
(warn "Array out of bonds for tridiagonal matrix. Ignored.")))))
(defmethod (setf vref) (value (matrix matrix-base-dgt) idx)
@@ -121,6 +120,6 @@
((< idx (- (* 3 len) 2))
(setf (aref (slot-value matrix 'subdiagonal-store) (- idx (- (* 2 len) 1)))
val2))
- (t
+ #+nil (t
(warn "Array out of bonds for tridiagonal matrix. Ignored.")))
val2))
More information about the lisplab-cvs
mailing list