[lisplab-cvs] r127 - src/matlisp
Jørn Inge Vestgården
jivestgarden at common-lisp.net
Thu Jan 14 14:14:32 UTC 2010
Author: jivestgarden
Date: Thu Jan 14 09:14:30 2010
New Revision: 127
Log:
contribution with potrf and potrs. Callable but not verified
Added:
src/matlisp/potrf.lisp
src/matlisp/potrs.lisp
Modified:
src/matlisp/lapack.lisp
Modified: src/matlisp/lapack.lisp
==============================================================================
--- src/matlisp/lapack.lisp (original)
+++ src/matlisp/lapack.lisp Thu Jan 14 09:14:30 2010
@@ -1866,5 +1866,281 @@
(LDB :integer :input)
(INFO :integer :output))
+(def-fortran-routine dpotrf :void
+ "
+ SUBROUTINE DPOTRF( UPLO, N, A, LDA, INFO )
+
+ -- LAPACK routine (version 3.1) --
+ Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
+ November 2006
+
+ .. Scalar Arguments ..
+ CHARACTER UPLO
+ INTEGER INFO, LDA, N
+ ..
+ .. Array Arguments ..
+ DOUBLE PRECISION A( LDA, * )
+ ..
+
+ Purpose
+ =======
+
+ DPOTRF computes the Cholesky factorization of a real symmetric
+ positive definite matrix A.
+
+ The factorization has the form
+ A = U**T * U, if UPLO = 'U', or
+ A = L * L**T, if UPLO = 'L',
+ where U is an upper triangular matrix and L is lower triangular.
+
+ This is the block version of the algorithm, calling Level 3 BLAS.
+
+ Arguments
+ =========
+
+ UPLO (input) CHARACTER*1
+ = 'U': Upper triangle of A is stored;
+ = 'L': Lower triangle of A is stored.
+
+ N (input) INTEGER
+ The order of the matrix A. N >= 0.
+
+ A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
+ On entry, the symmetric matrix A. If UPLO = 'U', the leading
+ N-by-N upper triangular part of A contains the upper
+ triangular part of the matrix A, and the strictly lower
+ triangular part of A is not referenced. If UPLO = 'L', the
+ leading N-by-N lower triangular part of A contains the lower
+ triangular part of the matrix A, and the strictly upper
+ triangular part of A is not referenced.
+
+ On exit, if INFO = 0, the factor U or L from the Cholesky
+ factorization A = U**T*U or A = L*L**T.
+
+ LDA (input) INTEGER
+ The leading dimension of the array A. LDA >= 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, the leading minor of order i is not
+ positive definite, and the factorization could not be
+ completed.
+
+ =====================================================================
+
+"
+ (uplo :string :input)
+ (n :integer :input)
+ (a (* :double-float) :input-output)
+ (lda :integer :input)
+ (info :integer :output))
+
+
+(def-fortran-routine dpotrs :void
+ "
+ SUBROUTINE DPOTRS( UPLO, N, NRHS, A, LDA, B, LDB, INFO )
+
+ -- LAPACK routine (version 3.1) --
+ Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
+ November 2006
+
+ .. Scalar Arguments ..
+ CHARACTER UPLO
+ INTEGER INFO, LDA, LDB, N, NRHS
+ ..
+ .. Array Arguments ..
+ DOUBLE PRECISION A( LDA, * ), B( LDB, * )
+ ..
+
+ Purpose
+ =======
+
+ DPOTRS solves a system of linear equations A*X = B with a symmetric
+ positive definite matrix A using the Cholesky factorization
+ A = U**T*U or A = L*L**T computed by DPOTRF.
+
+ Arguments
+ =========
+
+ UPLO (input) CHARACTER*1
+ = 'U': Upper triangle of A is stored;
+ = 'L': Lower triangle of A is stored.
+
+ 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 triangular factor U or L from the Cholesky factorization
+ A = U**T*U or A = L*L**T, as computed by DPOTRF.
+
+ LDA (input) INTEGER
+ The leading dimension of the array A. LDA >= max(1,N).
+
+ 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
+
+ =====================================================================
+
+"
+ (uplo :string :input)
+ (n :integer :input)
+ (nrhs :integer :input)
+ (a (* :double-float) :input)
+ (lda :integer :input)
+ (b (* :double-float) :input-output)
+ (ldb :integer :input)
+ (info :integer :output))
+
+
+(def-fortran-routine zpotrf :void
+ "
+ SUBROUTINE ZPOTRF( UPLO, N, A, LDA, INFO )
+
+ -- LAPACK routine (version 3.1) --
+ Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
+ November 2006
+
+ .. Scalar Arguments ..
+ CHARACTER UPLO
+ INTEGER INFO, LDA, N
+ ..
+ .. Array Arguments ..
+ COMPLEX*16 A( LDA, * )
+ ..
+
+ Purpose
+ =======
+
+ ZPOTRF computes the Cholesky factorization of a complex Hermitian
+ positive definite matrix A.
+
+ The factorization has the form
+ A = U**H * U, if UPLO = 'U', or
+ A = L * L**H, if UPLO = 'L',
+ where U is an upper triangular matrix and L is lower triangular.
+
+ This is the block version of the algorithm, calling Level 3 BLAS.
+
+ Arguments
+ =========
+
+ UPLO (input) CHARACTER*1
+ = 'U': Upper triangle of A is stored;
+ = 'L': Lower triangle of A is stored.
+
+ N (input) INTEGER
+ The order of the matrix A. N >= 0.
+
+ A (input/output) COMPLEX*16 array, dimension (LDA,N)
+ On entry, the Hermitian matrix A. If UPLO = 'U', the leading
+ N-by-N upper triangular part of A contains the upper
+ triangular part of the matrix A, and the strictly lower
+ triangular part of A is not referenced. If UPLO = 'L', the
+ leading N-by-N lower triangular part of A contains the lower
+ triangular part of the matrix A, and the strictly upper
+ triangular part of A is not referenced.
+
+ On exit, if INFO = 0, the factor U or L from the Cholesky
+ factorization A = U**H*U or A = L*L**H.
+
+ LDA (input) INTEGER
+ The leading dimension of the array A. LDA >= 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, the leading minor of order i is not
+ positive definite, and the factorization could not be
+ completed.
+
+ =====================================================================
+
+"
+ (uplo :string :input)
+ (n :integer :input)
+ (a (* :complex-double-float) :input-output)
+ (lda :integer :input)
+ (info :integer :output))
+
+
+(def-fortran-routine zpotrs :void
+ "
+ SUBROUTINE ZPOTRS( UPLO, N, NRHS, A, LDA, B, LDB, INFO )
+
+ -- LAPACK routine (version 3.1) --
+ Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
+ November 2006
+
+ .. Scalar Arguments ..
+ CHARACTER UPLO
+ INTEGER INFO, LDA, LDB, N, NRHS
+ ..
+ .. Array Arguments ..
+ COMPLEX*16 A( LDA, * ), B( LDB, * )
+ ..
+
+ Purpose
+ =======
+
+ ZPOTRS solves a system of linear equations A*X = B with a Hermitian
+ positive definite matrix A using the Cholesky factorization
+ A = U**H*U or A = L*L**H computed by ZPOTRF.
+
+ Arguments
+ =========
+
+ UPLO (input) CHARACTER*1
+ = 'U': Upper triangle of A is stored;
+ = 'L': Lower triangle of A is stored.
+
+ 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 triangular factor U or L from the Cholesky factorization
+ A = U**H*U or A = L*L**H, as computed by ZPOTRF.
+
+ LDA (input) INTEGER
+ The leading dimension of the array A. LDA >= max(1,N).
+
+ 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
+
+ =====================================================================
+
+"
+ (uplo :string :input)
+ (n :integer :input)
+ (nrhs :integer :input)
+ (a (* :complex-double-float) :input)
+ (lda :integer :input)
+ (b (* :complex-double-float) :input-output)
+ (ldb :integer :input)
+ (info :integer :output))
Added: src/matlisp/potrf.lisp
==============================================================================
--- (empty file)
+++ src/matlisp/potrf.lisp Thu Jan 14 09:14:30 2010
@@ -0,0 +1,92 @@
+;;; -*- Mode: lisp; Syntax: ansi-common-lisp; Package: :lisplab; Base: 10 -*-
+
+;;; Copyright (C) 2009 Knut S. Gjerden
+;;;
+;;; 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.
+
+;;; $Log: potrf.lisp,v $
+;;; Revision 1.2 14.01.2010 13:28:40 knutgj
+;;; Revision 1.1 06.08.2009 12:40:40 knutgj
+;;; o Initial revision.
+
+
+(in-package :lisplab)
+
+(defgeneric potrf! (a &key uplo)
+ (:documentation
+"
+ Syntax
+ ======
+ (POTRF a [:u :L])
+
+ Purpose
+ =======
+ DPOTRF computes the Cholesky factorization of a real symmetric
+ positive definite matrix A.
+
+ The factorization has the form
+ A = U**T * U, if UPLO = 'U', or
+ A = L * L**T, if UPLO = 'L',
+ where U is an upper triangular matrix and L is lower triangular.
+
+ This is the block version of the algorithm, calling Level 3 BLAS.
+
+ Return Values
+ =============
+ [1] The factor U or L from the Cholesky
+ factorization A = U**T*U or A = L*L**T.
+ [2] INFO = T: successful
+ i: U(i,i) is exactly zero.
+"))
+
+
+
+(defmethod potrf! ((a matrix-blas-dge) &key uplo)
+ (let* ((n (rows a))
+ (m (cols a)))
+
+ (declare (type fixnum n m))
+ (multiple-value-bind (new-a info)
+ (f77::dpotrf (case uplo
+ (:L "L")
+ (:U "U")
+ (t "U")) ;; UPLO
+ m ;; N
+ (matrix-store a) ;; A
+ n ;; LDA
+ 0) ;; INFO
+ (declare (ignore new-a))
+ (values a (if (zerop info)
+ t
+ info)))))
+
+(defmethod potrf! ((a matrix-blas-zge) &key uplo)
+ (let* ((n (rows a))
+ (m (cols a)))
+
+ (declare (type fixnum n m))
+ (multiple-value-bind (new-a info)
+ (f77::zpotrf (case uplo
+ (:L "L")
+ (:U "U")
+ (t "U")) ;; UPLO
+ m ;; N
+ (matrix-store a) ;; A
+ n ;; LDA
+ 0) ;; INFO
+ (declare (ignore new-a))
+ (values a (if (zerop info)
+ t
+ info)))))
Added: src/matlisp/potrs.lisp
==============================================================================
--- (empty file)
+++ src/matlisp/potrs.lisp Thu Jan 14 09:14:30 2010
@@ -0,0 +1,121 @@
+;;; -*- Mode: lisp; Syntax: ansi-common-lisp; Package: :lisplab; Base: 10 -*-
+
+;;; Copyright (C) 2009 Knut S. Gjerden
+;;;
+;;; 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.
+
+;;; $Log: potrs.lisp,v $
+;;; Revision 1.2 14.01.2010 13:29:27 knutgj
+;;; Revision 1.1 06.08.2009 12:19:27 knutgj
+;;; o Initial revision.
+
+
+(in-package :lisplab)
+
+(defgeneric potrs! (a b &key uplo)
+ (:documentation
+ "
+ Syntax
+ ======
+ (POTRS! a b [:U :L])
+
+ Purpose
+ =======
+ Solves a system of linear equations
+ A * X = B or A' * X = B
+ with a general N-by-N matrix A using the Cholesky LU factorization computed
+ by POTRF. A and are the results from POTRF, UPLO specifies
+ the form of the system of equations:
+ = 'U': A = U**T*U
+ = 'L': A = L*L**T
+
+ 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 potrs (a b &key uplo)
+ (:documentation
+ "
+ Syntax
+ ======
+ (POTRS a b [:U :L])
+
+ Purpose
+ =======
+ Same as POTRS! except that B is not overwritten.
+"))
+
+(defmethod potrs! :before ((a matrix-blas-dge) (b matrix-blas-dge) &key uplo)
+ (declare (ignore uplo))
+ (let ((n-a (rows a))
+ (m-a (cols a))
+ (n-b (rows b)))
+ (if (not (= n-a m-a n-b))
+ (error "Dimensions of A,B given to POTRS! do not match"))))
+
+(defmethod potrs! ((a matrix-blas-dge) (b matrix-blas-dge) &key uplo)
+ (let* ((n (rows a))
+ (m (cols b)))
+ (multiple-value-bind (x info)
+ (f77::dpotrs (case uplo
+ (:L "L")
+ (:U "U")
+ (t "U")) ;;UPLO
+ n ;;N
+ m ;;NRHS
+ (matrix-store a) ;;A
+ n ;;LDA
+ (matrix-store b) ;;B
+ n ;;LDB
+ 0) ;;INFO
+ (values
+ (make-instance 'matrix-dge :dim (list n m) :store x)
+ (if (zerop info)
+ t
+ info)))))
+
+(defmethod potrs! ((a matrix-blas-zge) (b matrix-blas-zge) &key uplo)
+
+ (let* ((n (rows a))
+ (m (cols b)))
+ (multiple-value-bind (x info)
+ (f77::zpotrs (case uplo
+ (:L "L")
+ (:U "U")
+ (t "U"))
+ n
+ m
+ (matrix-store a)
+ n
+ (matrix-store b)
+ n
+ 0)
+ (values
+ (make-instance 'matrix-zge :dim (list n m) :store x)
+ (if (zerop info)
+ t
+ info)))))
+
+(defmethod potrs ((a matrix-blas-dge) (b matrix-blas-dge) &key uplo)
+ (potrs! a (copy b) :uplo uplo))
+
+(defmethod potrs ((a matrix-blas-zge) (b matrix-blas-zge) &key uplo)
+ (potrs! a (copy b) :uplo uplo))
More information about the lisplab-cvs
mailing list