/*
* Copyright (c) 2002-2006 Samit Basu
*
* 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., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*
*/
#include "LeastSquaresSolver.hpp"
#include "LAPACK.hpp"
#include <stdlib.h>
#include <stdio.h>
#include "Malloc.hpp"
#define MSGBUFLEN 2048
/***************************************************************************
* Least-squares solver for double matrices
***************************************************************************/
/**
* Solve A * X = B in a least-squares sense, where A is m x n, and B is m x k.
* C is n x k.
*/
void doubleSolveLeastSq(Interpreter* eval,int m, int n, int k, double *c,
double *a, double *b) {
if ((m == 0) || (n == 0)) return;
char msgBuffer[MSGBUFLEN];
// Here are the comments from the LAPACK routine used:
//SUBROUTINE DGELSY( M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK,
// WORK, LWORK, INFO )
//*
//* -- LAPACK driver routine (version 3.0) --
//* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
//* Courant Institute, Argonne National Lab, and Rice University
//* June 30, 1999
//*
//* .. Scalar Arguments ..
// INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS, RANK
// DOUBLE PRECISION RCOND
//* ..
//* .. Array Arguments ..
// INTEGER JPVT( * )
// DOUBLE PRECISION A( LDA, * ), B( LDB, * ), WORK( * )
//* ..
//*
//* Purpose
//* =======
//*
//* DGELSY computes the minimum-norm solution to a real linear least
//* squares problem:
//* minimize || A * X - B ||
//* using a complete orthogonal factorization of A. A is an M-by-N
//* matrix which may be rank-deficient.
//*
//* Several right hand side vectors b and solution vectors x can be
//* handled in a single call; they are stored as the columns of the
//* M-by-NRHS right hand side matrix B and the N-by-NRHS solution
//* matrix X.
//*
//* The routine first computes a QR factorization with column pivoting:
//* A * P = Q * [ R11 R12 ]
//* [ 0 R22 ]
//* with R11 defined as the largest leading submatrix whose estimated
//* condition number is less than 1/RCOND. The order of R11, RANK,
//* is the effective rank of A.
//*
//* Then, R22 is considered to be negligible, and R12 is annihilated
//* by orthogonal transformations from the right, arriving at the
//* complete orthogonal factorization:
//* A * P = Q * [ T11 0 ] * Z
//* [ 0 0 ]
//* The minimum-norm solution is then
//* X = P * Z' [ inv(T11)*Q1'*B ]
//* [ 0 ]
//* where Q1 consists of the first RANK columns of Q.
//*
//* This routine is basically identical to the original xGELSX except
//* three differences:
//* o The call to the subroutine xGEQPF has been substituted by the
//* the call to the subroutine xGEQP3. This subroutine is a Blas-3
//* version of the QR factorization with column pivoting.
//* o Matrix B (the right hand side) is updated with Blas-3.
//* o The permutation of matrix B (the right hand side) is faster and
//* more simple.
//*
//* Arguments
//* =========
//*
//* M (input) INTEGER
//* The number of rows of the matrix A. M >= 0.
int M = m;
//* N (input) INTEGER
//* The number of columns of the matrix A. N >= 0.
int N = n;
//* NRHS (input) INTEGER
//* The number of right hand sides, i.e., the number of
//* columns of matrices B and X. NRHS >= 0.
int NRHS = k;
//* A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
//* On entry, the M-by-N matrix A.
//* On exit, A has been overwritten by details of its
//* complete orthogonal factorization.
double *A = a;
//* LDA (input) INTEGER
//* The leading dimension of the array A. LDA >= max(1,M).
int LDA = m;
//* B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
//* On entry, the M-by-NRHS right hand side matrix B.
//* On exit, the N-by-NRHS solution matrix X.
double *B;
int Bsize = (M > N) ? M : N;
// This passing convention requires that we copy our source matrix
// into the destination array with the appropriate padding.
B = (double*) Calloc(Bsize*NRHS*sizeof(double));
changeStrideDouble(B,Bsize,b,m,m,NRHS);
//* LDB (input) INTEGER
//* The leading dimension of the array B. LDB >= max(1,M,N).
int LDB = Bsize;
//* JPVT (input/output) INTEGER array, dimension (N)
//* On entry, if JPVT(i) .ne. 0, the i-th column of A is permuted
//* to the front of AP, otherwise column i is a Free column.
//* On exit, if JPVT(i) = k, then the i-th column of AP
//* was the k-th column of A.
int *JPVT = (int*) Calloc(N*sizeof(int));
//* RCOND (input) DOUBLE PRECISION
//* RCOND is used to determine the effective rank of A, which
//* is defined as the order of the largest leading triangular
//* submatrix R11 in the QR factorization with pivoting of A,
//* whose estimated condition number < 1/RCOND.
double RCOND = getEPS();
//* RANK (output) INTEGER
//* The effective rank of A, i.e., the order of the submatrix
//* R11. This is the same as the order of the submatrix T11
//* in the complete orthogonal factorization of A.
int RANK;
//* WORK (workspace/output) DOUBLE PRECISION array, dimension (LWORK)
//* On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
double WORKSIZE;
//* LWORK (input) INTEGER
//* The dimension of the array WORK.
//* The unblocked strategy requires that:
//* LWORK >= MAX( MN+3*N+1, 2*MN+NRHS ),
//* where MN = min( M, N ).
//* The block algorithm requires that:
//* LWORK >= MAX( MN+2*N+NB*(N+1), 2*MN+NB*NRHS ),
//* where NB is an upper bound on the blocksize returned
//* by ILAENV for the routines DGEQP3, DTZRZF, STZRQF, DORMQR,
//* and DORMRZ.
//*
//* If LWORK = -1, then a workspace query is assumed; the routine
//* only calculates the optimal size of the WORK array, returns
//* this value as the first entry of the WORK array, and no error
//* message related to LWORK is issued by XERBLA.
int LWORK;
//* INFO (output) INTEGER
//* = 0: successful exit
//* < 0: If INFO = -i, the i-th argument had an illegal value.
int INFO;
//* Further Details
//* ===============
//*
//* Based on contributions by
//* A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA
//* E. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain
//* G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain
//*
//* =====================================================================
//
LWORK = -1;
dgelsy_(&M, &N, &NRHS, A, &LDA, B, &LDB, JPVT, &RCOND,
&RANK, &WORKSIZE, &LWORK, &INFO);
LWORK = (int) WORKSIZE;
double *WORK = (double*) Malloc(LWORK*sizeof(double));
dgelsy_(&M, &N, &NRHS, A, &LDA, B, &LDB, JPVT, &RCOND,
&RANK, WORK, &LWORK, &INFO);
// Check the rank...
if (M > N) {
// Problem should be overdetermined, rank should be N
if (RANK < N) {
snprintf(msgBuffer,MSGBUFLEN,"Matrix is rank deficient to machine precision. RANK = %d\n",RANK);
eval->warningMessage(msgBuffer);
}
} else
// Problem should be underderemined, rank should be M
if (RANK < M) {
snprintf(msgBuffer,MSGBUFLEN,"Matrix is rank deficient to machine precision. RANK = %d\n",RANK);
eval->warningMessage(msgBuffer);
}
changeStrideDouble(c,n,B,Bsize,n,k);
// Free the allocated arrays
Free(B);
Free(JPVT);
Free(WORK);
}
/***************************************************************************
* Least-squares solver for complex matrices
***************************************************************************/
/**
* Solve A * X = B in a least-squares sense, where A is m x n, and B is m x k.
* C is n x k.
*/
void dcomplexSolveLeastSq(Interpreter* eval,int m, int n, int k, double *c,
double *a, double *b) {
if ((m == 0) || (n == 0)) return;
char msgBuffer[MSGBUFLEN];
// SUBROUTINE ZGELSY( M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK,
// $ WORK, LWORK, RWORK, INFO )
//*
//* -- LAPACK driver routine (version 3.0) --
//* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
//* Courant Institute, Argonne National Lab, and Rice University
//* June 30, 1999
//*
//* .. Scalar Arguments ..
// INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS, RANK
// DOUBLE PRECISION RCOND
//* ..
//* .. Array Arguments ..
// INTEGER JPVT( * )
// DOUBLE PRECISION RWORK( * )
// COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * )
//* ..
//*
//* Purpose
//* =======
//*
//* ZGELSY computes the minimum-norm solution to a real linear least
//* squares problem:
//* minimize || A * X - B ||
//* using a complete orthogonal factorization of A. A is an M-by-N
//* matrix which may be rank-deficient.
//*
//* Several right hand side vectors b and solution vectors x can be
//* handled in a single call; they are stored as the columns of the
//* M-by-NRHS right hand side matrix B and the N-by-NRHS solution
//* matrix X.
//*
//* The routine first computes a QR factorization with column pivoting:
//* A * P = Q * [ R11 R12 ]
//* [ 0 R22 ]
//* with R11 defined as the largest leading submatrix whose estimated
//* condition number is less than 1/RCOND. The order of R11, RANK,
//* is the effective rank of A.
//*
//* Then, R22 is considered to be negligible, and R12 is annihilated
//* by orthogonal transformations from the right, arriving at the
//* complete orthogonal factorization:
//* A * P = Q * [ T11 0 ] * Z
//* [ 0 0 ]
//* The minimum-norm solution is then
//* X = P * Z' [ inv(T11)*Q1'*B ]
//* [ 0 ]
//* where Q1 consists of the first RANK columns of Q.
//*
//* This routine is basically identical to the original xGELSX except
//* three differences:
//* o The call to the subroutine xGEQPF has been substituted by the
//* the call to the subroutine xGEQP3. This subroutine is a Blas-3
//* version of the QR factorization with column pivoting.
//* o Matrix B (the right hand side) is updated with Blas-3.
//* o The permutation of matrix B (the right hand side) is faster and
//* more simple.
//*
//* Arguments
//* =========
//*
//* M (input) INTEGER
//* The number of rows of the matrix A. M >= 0.
int M = m;
//* N (input) INTEGER
//* The number of columns of the matrix A. N >= 0.
int N = n;
//* NRHS (input) INTEGER
//* The number of right hand sides, i.e., the number of
//* columns of matrices B and X. NRHS >= 0.
int NRHS = k;
//* A (input/output) COMPLEX*16 array, dimension (LDA,N)
//* On entry, the M-by-N matrix A.
//* On exit, A has been overwritten by details of its
//* complete orthogonal factorization.
double *A = a;
//* LDA (input) INTEGER
//* The leading dimension of the array A. LDA >= max(1,M).
int LDA = m;
//* B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
//* On entry, the M-by-NRHS right hand side matrix B.
//* On exit, the N-by-NRHS solution matrix X.
double *B;
int Bsize = (M > N) ? M : N;
// This passing convention requires that we copy our source matrix
// into the destination array with the appropriate padding.
B = (double*) Calloc(Bsize*NRHS*2*sizeof(double));
changeStrideDouble(B,2*Bsize,b,2*m,2*m,NRHS);
//* LDB (input) INTEGER
//* The leading dimension of the array B. LDB >= max(1,M,N).
int LDB = Bsize;
//* JPVT (input/output) INTEGER array, dimension (N)
//* On entry, if JPVT(i) .ne. 0, the i-th column of A is permuted
//* to the front of AP, otherwise column i is a Free column.
//* On exit, if JPVT(i) = k, then the i-th column of AP
//* was the k-th column of A.
int *JPVT = (int*) Calloc(N*sizeof(int));
//* RCOND (input) DOUBLE PRECISION
//* RCOND is used to determine the effective rank of A, which
//* is defined as the order of the largest leading triangular
//* submatrix R11 in the QR factorization with pivoting of A,
//* whose estimated condition number < 1/RCOND.
double RCOND = getEPS();
//* RANK (output) INTEGER
//* The effective rank of A, i.e., the order of the submatrix
//* R11. This is the same as the order of the submatrix T11
//* in the complete orthogonal factorization of A.
int RANK;
//* WORK (workspace/output) COMPLEX*16 array, dimension (LWORK)
//* On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
double WORKSIZE;
//* LWORK (input) INTEGER
//* The dimension of the array WORK.
//* The unblocked strategy requires that:
//* LWORK >= MN + MAX( 2*MN, N+1, MN+NRHS )
//* where MN = min(M,N).
//* The block algorithm requires that:
//* LWORK >= MN + MAX( 2*MN, NB*(N+1), MN+MN*NB, MN+NB*NRHS )
//* where NB is an upper bound on the blocksize returned
//* by ILAENV for the routines ZGEQP3, ZTZRZF, CTZRQF, ZUNMQR,
//* and ZUNMRZ.
//*
//* If LWORK = -1, then a workspace query is assumed; the routine
//* only calculates the optimal size of the WORK array, returns
//* this value as the first entry of the WORK array, and no error
//* message related to LWORK is issued by XERBLA.
int LWORK;
//* RWORK (workspace) DOUBLE PRECISION array, dimension (2*N)
double *RWORK = (double*) Malloc(2*N*sizeof(double));
//* INFO (output) INTEGER
//* = 0: successful exit
//* < 0: If INFO = -i, the i-th argument had an illegal value.
int INFO;
//* Further Details
//* ===============
//*
//* Based on contributions by
//* A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA
//* E. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain
//* G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain
//*
//* =====================================================================
//
LWORK = -1;
zgelsy_(&M, &N, &NRHS, A, &LDA, B, &LDB, JPVT, &RCOND,
&RANK, &WORKSIZE, &LWORK, RWORK, &INFO);
LWORK = (int) WORKSIZE;
double *WORK = (double*) Malloc(LWORK*sizeof(double));
zgelsy_(&M, &N, &NRHS, A, &LDA, B, &LDB, JPVT, &RCOND,
&RANK, WORK, &LWORK, RWORK, &INFO);
// Check the rank...
if (M > N) {
// Problem should be overdetermined, rank should be N
if (RANK < N)
snprintf(msgBuffer,MSGBUFLEN,"Matrix is rank deficient to machine precision. RANK = %d\n",RANK);
eval->warningMessage(msgBuffer);
} else
// Problem should be underderemined, rank should be M
if (RANK < M) {
snprintf(msgBuffer,MSGBUFLEN,"Matrix is rank deficient to machine precision. RANK = %d\n",RANK);
eval->warningMessage(msgBuffer);
}
changeStrideDouble(c,2*n,B,2*Bsize,2*n,k);
// Free the allocated arrays
Free(B);
Free(JPVT);
Free(WORK);
Free(RWORK);
}
/***************************************************************************
* Least-squares solver for float matrices
***************************************************************************/
/**
* Solve A * X = B in a least-squares sense, where A is m x n, and B is m x k.
* C is n x k.
*/
void floatSolveLeastSq(Interpreter* eval,int m, int n, int k, float *c,
float *a, float *b) {
if ((m == 0) || (n == 0)) return;
char msgBuffer[MSGBUFLEN];
// Here are the comments from the LAPACK routine used:
//SUBROUTINE SGELSY( M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK,
// WORK, LWORK, INFO )
//*
//* -- LAPACK driver routine (version 3.0) --
//* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
//* Courant Institute, Argonne National Lab, and Rice University
//* June 30, 1999
//*
//* .. Scalar Arguments ..
// INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS, RANK
// REAL RCOND
//* ..
//* .. Array Arguments ..
// INTEGER JPVT( * )
// REAL A( LDA, * ), B( LDB, * ), WORK( * )
//* ..
//*
//* Purpose
//* =======
//*
//* SGELSY computes the minimum-norm solution to a real linear least
//* squares problem:
//* minimize || A * X - B ||
//* using a complete orthogonal factorization of A. A is an M-by-N
//* matrix which may be rank-deficient.
//*
//* Several right hand side vectors b and solution vectors x can be
//* handled in a single call; they are stored as the columns of the
//* M-by-NRHS right hand side matrix B and the N-by-NRHS solution
//* matrix X.
//*
//* The routine first computes a QR factorization with column pivoting:
//* A * P = Q * [ R11 R12 ]
//* [ 0 R22 ]
//* with R11 defined as the largest leading submatrix whose estimated
//* condition number is less than 1/RCOND. The order of R11, RANK,
//* is the effective rank of A.
//*
//* Then, R22 is considered to be negligible, and R12 is annihilated
//* by orthogonal transformations from the right, arriving at the
//* complete orthogonal factorization:
//* A * P = Q * [ T11 0 ] * Z
//* [ 0 0 ]
//* The minimum-norm solution is then
//* X = P * Z' [ inv(T11)*Q1'*B ]
//* [ 0 ]
//* where Q1 consists of the first RANK columns of Q.
//*
//* This routine is basically identical to the original xGELSX except
//* three differences:
//* o The call to the subroutine xGEQPF has been substituted by the
//* the call to the subroutine xGEQP3. This subroutine is a Blas-3
//* version of the QR factorization with column pivoting.
//* o Matrix B (the right hand side) is updated with Blas-3.
//* o The permutation of matrix B (the right hand side) is faster and
//* more simple.
//*
//* Arguments
//* =========
//*
//* M (input) INTEGER
//* The number of rows of the matrix A. M >= 0.
int M = m;
//* N (input) INTEGER
//* The number of columns of the matrix A. N >= 0.
int N = n;
//* NRHS (input) INTEGER
//* The number of right hand sides, i.e., the number of
//* columns of matrices B and X. NRHS >= 0.
int NRHS = k;
//* A (input/output) REAL array, dimension (LDA,N)
//* On entry, the M-by-N matrix A.
//* On exit, A has been overwritten by details of its
//* complete orthogonal factorization.
float *A = a;
//* LDA (input) INTEGER
//* The leading dimension of the array A. LDA >= max(1,M).
int LDA = m;
//* B (input/output) REAL array, dimension (LDB,NRHS)
//* On entry, the M-by-NRHS right hand side matrix B.
//* On exit, the N-by-NRHS solution matrix X.
float *B;
int Bsize = (M > N) ? M : N;
// This passing convention requires that we copy our source matrix
// into the destination array with the appropriate padding.
B = (float*) Calloc(Bsize*NRHS*sizeof(float));
changeStrideFloat(B,Bsize,b,m,m,NRHS);
//* LDB (input) INTEGER
//* The leading dimension of the array B. LDB >= max(1,M,N).
int LDB = Bsize;
//* JPVT (input/output) INTEGER array, dimension (N)
//* On entry, if JPVT(i) .ne. 0, the i-th column of A is permuted
//* to the front of AP, otherwise column i is a Free column.
//* On exit, if JPVT(i) = k, then the i-th column of AP
//* was the k-th column of A.
int *JPVT = (int*) Calloc(N*sizeof(int));
//* RCOND (input) REAL
//* RCOND is used to determine the effective rank of A, which
//* is defined as the order of the largest leading triangular
//* submatrix R11 in the QR factorization with pivoting of A,
//* whose estimated condition number < 1/RCOND.
float RCOND = getFloatEPS();
//* RANK (output) INTEGER
//* The effective rank of A, i.e., the order of the submatrix
//* R11. This is the same as the order of the submatrix T11
//* in the complete orthogonal factorization of A.
int RANK;
//* WORK (workspace/output) REAL array, dimension (LWORK)
//* On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
float WORKSIZE;
//* LWORK (input) INTEGER
//* The dimension of the array WORK.
//* The unblocked strategy requires that:
//* LWORK >= MAX( MN+3*N+1, 2*MN+NRHS ),
//* where MN = min( M, N ).
//* The block algorithm requires that:
//* LWORK >= MAX( MN+2*N+NB*(N+1), 2*MN+NB*NRHS ),
//* where NB is an upper bound on the blocksize returned
//* by ILAENV for the routines DGEQP3, DTZRZF, STZRQF, DORMQR,
//* and DORMRZ.
//*
//* If LWORK = -1, then a workspace query is assumed; the routine
//* only calculates the optimal size of the WORK array, returns
//* this value as the first entry of the WORK array, and no error
//* message related to LWORK is issued by XERBLA.
int LWORK;
//* INFO (output) INTEGER
//* = 0: successful exit
//* < 0: If INFO = -i, the i-th argument had an illegal value.
int INFO;
//* Further Details
//* ===============
//*
//* Based on contributions by
//* A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA
//* E. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain
//* G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain
//*
//* =====================================================================
//
LWORK = -1;
sgelsy_(&M, &N, &NRHS, A, &LDA, B, &LDB, JPVT, &RCOND,
&RANK, &WORKSIZE, &LWORK, &INFO);
LWORK = (int) WORKSIZE;
float *WORK = (float*) Malloc(LWORK*sizeof(float));
sgelsy_(&M, &N, &NRHS, A, &LDA, B, &LDB, JPVT, &RCOND,
&RANK, WORK, &LWORK, &INFO);
// Check the rank...
if (M > N) {
// Problem should be overdetermined, rank should be N
if (RANK < N) {
snprintf(msgBuffer,MSGBUFLEN,"Matrix is rank deficient to machine precision. RANK = %d\n",RANK);
eval->warningMessage(msgBuffer);
}
} else
// Problem should be underderemined, rank should be M
if (RANK < M) {
snprintf(msgBuffer,MSGBUFLEN,"Matrix is rank deficient to machine precision. RANK = %d\n",RANK);
eval->warningMessage(msgBuffer);
}
changeStrideFloat(c,n,B,Bsize,n,k);
// Free the allocated arrays
Free(B);
Free(JPVT);
Free(WORK);
}
/***************************************************************************
* Least-squares solver for complex matrices
***************************************************************************/
/**
* Solve A * X = B in a least-squares sense, where A is m x n, and B is m x k.
* C is n x k.
*/
void complexSolveLeastSq(Interpreter* eval,int m, int n, int k, float *c,
float *a, float *b) {
if ((m == 0) || (n == 0)) return;
char msgBuffer[MSGBUFLEN];
// SUBROUTINE CGELSY( M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK,
// $ WORK, LWORK, RWORK, INFO )
//*
//* -- LAPACK driver routine (version 3.0) --
//* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
//* Courant Institute, Argonne National Lab, and Rice University
//* June 30, 1999
//*
//* .. Scalar Arguments ..
// INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS, RANK
// REAL RCOND
//* ..
//* .. Array Arguments ..
// INTEGER JPVT( * )
// REAL RWORK( * )
// COMPLEX A( LDA, * ), B( LDB, * ), WORK( * )
//* ..
//*
//* Purpose
//* =======
//*
//* CGELSY computes the minimum-norm solution to a real linear least
//* squares problem:
//* minimize || A * X - B ||
//* using a complete orthogonal factorization of A. A is an M-by-N
//* matrix which may be rank-deficient.
//*
//* Several right hand side vectors b and solution vectors x can be
//* handled in a single call; they are stored as the columns of the
//* M-by-NRHS right hand side matrix B and the N-by-NRHS solution
//* matrix X.
//*
//* The routine first computes a QR factorization with column pivoting:
//* A * P = Q * [ R11 R12 ]
//* [ 0 R22 ]
//* with R11 defined as the largest leading submatrix whose estimated
//* condition number is less than 1/RCOND. The order of R11, RANK,
//* is the effective rank of A.
//*
//* Then, R22 is considered to be negligible, and R12 is annihilated
//* by orthogonal transformations from the right, arriving at the
//* complete orthogonal factorization:
//* A * P = Q * [ T11 0 ] * Z
//* [ 0 0 ]
//* The minimum-norm solution is then
//* X = P * Z' [ inv(T11)*Q1'*B ]
//* [ 0 ]
//* where Q1 consists of the first RANK columns of Q.
//*
//* This routine is basically identical to the original xGELSX except
//* three differences:
//* o The call to the subroutine xGEQPF has been substituted by the
//* the call to the subroutine xGEQP3. This subroutine is a Blas-3
//* version of the QR factorization with column pivoting.
//* o Matrix B (the right hand side) is updated with Blas-3.
//* o The permutation of matrix B (the right hand side) is faster and
//* more simple.
//*
//* Arguments
//* =========
//*
//* M (input) INTEGER
//* The number of rows of the matrix A. M >= 0.
int M = m;
//* N (input) INTEGER
//* The number of columns of the matrix A. N >= 0.
int N = n;
//* NRHS (input) INTEGER
//* The number of right hand sides, i.e., the number of
//* columns of matrices B and X. NRHS >= 0.
int NRHS = k;
//* A (input/output) COMPLEX array, dimension (LDA,N)
//* On entry, the M-by-N matrix A.
//* On exit, A has been overwritten by details of its
//* complete orthogonal factorization.
float *A = a;
//* LDA (input) INTEGER
//* The leading dimension of the array A. LDA >= max(1,M).
int LDA = m;
//* B (input/output) REAL array, dimension (LDB,NRHS)
//* On entry, the M-by-NRHS right hand side matrix B.
//* On exit, the N-by-NRHS solution matrix X.
float *B;
int Bsize = (M > N) ? M : N;
// This passing convention requires that we copy our source matrix
// into the destination array with the appropriate padding.
B = (float*) Calloc(Bsize*NRHS*2*sizeof(float));
changeStrideFloat(B,2*Bsize,b,2*m,2*m,NRHS);
//* LDB (input) INTEGER
//* The leading dimension of the array B. LDB >= max(1,M,N).
int LDB = Bsize;
//* JPVT (input/output) INTEGER array, dimension (N)
//* On entry, if JPVT(i) .ne. 0, the i-th column of A is permuted
//* to the front of AP, otherwise column i is a Free column.
//* On exit, if JPVT(i) = k, then the i-th column of AP
//* was the k-th column of A.
int *JPVT = (int*) Calloc(N*sizeof(int));
//* RCOND (input) REAL
//* RCOND is used to determine the effective rank of A, which
//* is defined as the order of the largest leading triangular
//* submatrix R11 in the QR factorization with pivoting of A,
//* whose estimated condition number < 1/RCOND.
float RCOND = getFloatEPS();
//* RANK (output) INTEGER
//* The effective rank of A, i.e., the order of the submatrix
//* R11. This is the same as the order of the submatrix T11
//* in the complete orthogonal factorization of A.
int RANK;
//* WORK (workspace/output) COMPLEX array, dimension (LWORK)
//* On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
float WORKSIZE;
//* LWORK (input) INTEGER
//* The dimension of the array WORK.
//* The unblocked strategy requires that:
//* LWORK >= MN + MAX( 2*MN, N+1, MN+NRHS )
//* where MN = min(M,N).
//* The block algorithm requires that:
//* LWORK >= MN + MAX( 2*MN, NB*(N+1), MN+MN*NB, MN+NB*NRHS )
//* where NB is an upper bound on the blocksize returned
//* by ILAENV for the routines ZGEQP3, ZTZRZF, CTZRQF, ZUNMQR,
//* and ZUNMRZ.
//*
//* If LWORK = -1, then a workspace query is assumed; the routine
//* only calculates the optimal size of the WORK array, returns
//* this value as the first entry of the WORK array, and no error
//* message related to LWORK is issued by XERBLA.
int LWORK;
//* RWORK (workspace) REAL array, dimension (2*N)
float *RWORK = (float*) Malloc(2*N*sizeof(float));
//* INFO (output) INTEGER
//* = 0: successful exit
//* < 0: If INFO = -i, the i-th argument had an illegal value.
int INFO;
//* Further Details
//* ===============
//*
//* Based on contributions by
//* A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA
//* E. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain
//* G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain
//*
//* =====================================================================
//
LWORK = -1;
cgelsy_(&M, &N, &NRHS, A, &LDA, B, &LDB, JPVT, &RCOND,
&RANK, &WORKSIZE, &LWORK, RWORK, &INFO);
LWORK = (int) WORKSIZE;
float *WORK = (float*) Malloc(LWORK*sizeof(float));
cgelsy_(&M, &N, &NRHS, A, &LDA, B, &LDB, JPVT, &RCOND,
&RANK, WORK, &LWORK, RWORK, &INFO);
// Check the rank...
if (M > N) {
// Problem should be overdetermined, rank should be N
if (RANK < N) {
snprintf(msgBuffer,MSGBUFLEN,"Matrix is rank deficient to machine precision. RANK = %d\n",RANK);
eval->warningMessage(msgBuffer);
}
} else
// Problem should be underderemined, rank should be M
if (RANK < M) {
snprintf(msgBuffer,MSGBUFLEN,"Matrix is rank deficient to machine precision. RANK = %d\n",RANK);
eval->warningMessage(msgBuffer);
}
changeStrideFloat(c,2*n,B,2*Bsize,2*n,k);
// Free the allocated arrays
Free(B);
Free(JPVT);
Free(WORK);
Free(RWORK);
}
syntax highlighted by Code2HTML, v. 0.9.1