|
Packit |
67cb25 |
/* eigen/schur.c
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* Copyright (C) 2006, 2007 Patrick Alken
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* This program is free software; you can redistribute it and/or modify
|
|
Packit |
67cb25 |
* it under the terms of the GNU General Public License as published by
|
|
Packit |
67cb25 |
* the Free Software Foundation; either version 3 of the License, or (at
|
|
Packit |
67cb25 |
* your option) any later version.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* This program is distributed in the hope that it will be useful, but
|
|
Packit |
67cb25 |
* WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
Packit |
67cb25 |
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
|
Packit |
67cb25 |
* General Public License for more details.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* You should have received a copy of the GNU General Public License
|
|
Packit |
67cb25 |
* along with this program; if not, write to the Free Software
|
|
Packit |
67cb25 |
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#include <config.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_eigen.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_math.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_matrix.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_vector.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_vector_complex.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_blas.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_complex.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_complex_math.h>
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
* This module contains some routines related to manipulating the
|
|
Packit |
67cb25 |
* Schur form of a matrix which are needed by the eigenvalue solvers
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* This file contains routines based on original code from LAPACK
|
|
Packit |
67cb25 |
* which is distributed under the modified BSD license.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#define GSL_SCHUR_SMLNUM (2.0 * GSL_DBL_MIN)
|
|
Packit |
67cb25 |
#define GSL_SCHUR_BIGNUM ((1.0 - GSL_DBL_EPSILON) / GSL_SCHUR_SMLNUM)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
gsl_schur_gen_eigvals()
|
|
Packit |
67cb25 |
Compute the eigenvalues of a 2-by-2 generalized block.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Inputs: A - 2-by-2 matrix
|
|
Packit |
67cb25 |
B - 2-by-2 upper triangular matrix
|
|
Packit |
67cb25 |
wr1 - (output) see notes
|
|
Packit |
67cb25 |
wr2 - (output) see notes
|
|
Packit |
67cb25 |
wi - (output) see notes
|
|
Packit |
67cb25 |
scale1 - (output) see notes
|
|
Packit |
67cb25 |
scale2 - (output) see notes
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Return: success
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Notes:
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
1)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
If the block contains real eigenvalues, then wi is set to 0,
|
|
Packit |
67cb25 |
and wr1, wr2, scale1, and scale2 are set such that:
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
eval1 = wr1 * scale1
|
|
Packit |
67cb25 |
eval2 = wr2 * scale2
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
If the block contains complex eigenvalues, then wr1, wr2, scale1,
|
|
Packit |
67cb25 |
scale2, and wi are set such that:
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
wr1 = wr2 = scale1 * Re(eval)
|
|
Packit |
67cb25 |
wi = scale1 * Im(eval)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
wi is always non-negative
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
2) This routine is based on LAPACK DLAG2
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_schur_gen_eigvals(const gsl_matrix *A, const gsl_matrix *B, double *wr1,
|
|
Packit |
67cb25 |
double *wr2, double *wi, double *scale1,
|
|
Packit |
67cb25 |
double *scale2)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const double safemin = GSL_DBL_MIN * 1.0e2;
|
|
Packit |
67cb25 |
const double safemax = 1.0 / safemin;
|
|
Packit |
67cb25 |
const double rtmin = sqrt(safemin);
|
|
Packit |
67cb25 |
const double rtmax = 1.0 / rtmin;
|
|
Packit |
67cb25 |
double anorm, bnorm;
|
|
Packit |
67cb25 |
double ascale, bscale, bsize;
|
|
Packit |
67cb25 |
double s1, s2;
|
|
Packit |
67cb25 |
double A11, A12, A21, A22;
|
|
Packit |
67cb25 |
double B11, B12, B22;
|
|
Packit |
67cb25 |
double binv11, binv22;
|
|
Packit |
67cb25 |
double bmin;
|
|
Packit |
67cb25 |
double as11, as12, as22, abi22;
|
|
Packit |
67cb25 |
double pp, qq, shift, ss, discr, r;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* scale A */
|
|
Packit |
67cb25 |
anorm = GSL_MAX(GSL_MAX(fabs(gsl_matrix_get(A, 0, 0)) +
|
|
Packit |
67cb25 |
fabs(gsl_matrix_get(A, 1, 0)),
|
|
Packit |
67cb25 |
fabs(gsl_matrix_get(A, 0, 1)) +
|
|
Packit |
67cb25 |
fabs(gsl_matrix_get(A, 1, 1))),
|
|
Packit |
67cb25 |
safemin);
|
|
Packit |
67cb25 |
ascale = 1.0 / anorm;
|
|
Packit |
67cb25 |
A11 = ascale * gsl_matrix_get(A, 0, 0);
|
|
Packit |
67cb25 |
A12 = ascale * gsl_matrix_get(A, 0, 1);
|
|
Packit |
67cb25 |
A21 = ascale * gsl_matrix_get(A, 1, 0);
|
|
Packit |
67cb25 |
A22 = ascale * gsl_matrix_get(A, 1, 1);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* perturb B if necessary to ensure non-singularity */
|
|
Packit |
67cb25 |
B11 = gsl_matrix_get(B, 0, 0);
|
|
Packit |
67cb25 |
B12 = gsl_matrix_get(B, 0, 1);
|
|
Packit |
67cb25 |
B22 = gsl_matrix_get(B, 1, 1);
|
|
Packit |
67cb25 |
bmin = rtmin * GSL_MAX(fabs(B11),
|
|
Packit |
67cb25 |
GSL_MAX(fabs(B12), GSL_MAX(fabs(B22), rtmin)));
|
|
Packit |
67cb25 |
if (fabs(B11) < bmin)
|
|
Packit |
67cb25 |
B11 = GSL_SIGN(B11) * bmin;
|
|
Packit |
67cb25 |
if (fabs(B22) < bmin)
|
|
Packit |
67cb25 |
B22 = GSL_SIGN(B22) * bmin;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* scale B */
|
|
Packit |
67cb25 |
bnorm = GSL_MAX(fabs(B11), GSL_MAX(fabs(B12) + fabs(B22), safemin));
|
|
Packit |
67cb25 |
bsize = GSL_MAX(fabs(B11), fabs(B22));
|
|
Packit |
67cb25 |
bscale = 1.0 / bsize;
|
|
Packit |
67cb25 |
B11 *= bscale;
|
|
Packit |
67cb25 |
B12 *= bscale;
|
|
Packit |
67cb25 |
B22 *= bscale;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* compute larger eigenvalue */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
binv11 = 1.0 / B11;
|
|
Packit |
67cb25 |
binv22 = 1.0 / B22;
|
|
Packit |
67cb25 |
s1 = A11 * binv11;
|
|
Packit |
67cb25 |
s2 = A22 * binv22;
|
|
Packit |
67cb25 |
if (fabs(s1) <= fabs(s2))
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
as12 = A12 - s1 * B12;
|
|
Packit |
67cb25 |
as22 = A22 - s1 * B22;
|
|
Packit |
67cb25 |
ss = A21 * (binv11 * binv22);
|
|
Packit |
67cb25 |
abi22 = as22 * binv22 - ss * B12;
|
|
Packit |
67cb25 |
pp = 0.5 * abi22;
|
|
Packit |
67cb25 |
shift = s1;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
as12 = A12 - s2 * B12;
|
|
Packit |
67cb25 |
as11 = A11 - s2 * B11;
|
|
Packit |
67cb25 |
ss = A21 * (binv11 * binv22);
|
|
Packit |
67cb25 |
abi22 = -ss * B12;
|
|
Packit |
67cb25 |
pp = 0.5 * (as11 * binv11 + abi22);
|
|
Packit |
67cb25 |
shift = s2;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
qq = ss * as12;
|
|
Packit |
67cb25 |
if (fabs(pp * rtmin) >= 1.0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
discr = (rtmin * pp) * (rtmin * pp) + qq * safemin;
|
|
Packit |
67cb25 |
r = sqrt(fabs(discr)) * rtmax;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (pp * pp + fabs(qq) <= safemin)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
discr = (rtmax * pp) * (rtmax * pp) + qq * safemax;
|
|
Packit |
67cb25 |
r = sqrt(fabs(discr)) * rtmin;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
discr = pp * pp + qq;
|
|
Packit |
67cb25 |
r = sqrt(fabs(discr));
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (discr >= 0.0 || r == 0.0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double sum = pp + GSL_SIGN(pp) * r;
|
|
Packit |
67cb25 |
double diff = pp - GSL_SIGN(pp) * r;
|
|
Packit |
67cb25 |
double wbig = shift + sum;
|
|
Packit |
67cb25 |
double wsmall = shift + diff;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* compute smaller eigenvalue */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (0.5 * fabs(wbig) > GSL_MAX(fabs(wsmall), safemin))
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double wdet = (A11*A22 - A12*A21) * (binv11 * binv22);
|
|
Packit |
67cb25 |
wsmall = wdet / wbig;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* choose (real) eigenvalue closest to 2,2 element of AB^{-1} for wr1 */
|
|
Packit |
67cb25 |
if (pp > abi22)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
*wr1 = GSL_MIN(wbig, wsmall);
|
|
Packit |
67cb25 |
*wr2 = GSL_MAX(wbig, wsmall);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
*wr1 = GSL_MAX(wbig, wsmall);
|
|
Packit |
67cb25 |
*wr2 = GSL_MIN(wbig, wsmall);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
*wi = 0.0;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* complex eigenvalues */
|
|
Packit |
67cb25 |
*wr1 = shift + pp;
|
|
Packit |
67cb25 |
*wr2 = *wr1;
|
|
Packit |
67cb25 |
*wi = r;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* compute scaling */
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const double fuzzy1 = 1.0 + 1.0e-5;
|
|
Packit |
67cb25 |
double c1, c2, c3, c4, c5;
|
|
Packit |
67cb25 |
double wabs, wsize, wscale;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
c1 = bsize * (safemin * GSL_MAX(1.0, ascale));
|
|
Packit |
67cb25 |
c2 = safemin * GSL_MAX(1.0, bnorm);
|
|
Packit |
67cb25 |
c3 = bsize * safemin;
|
|
Packit |
67cb25 |
if (ascale <= 1.0 && bsize <= 1.0)
|
|
Packit |
67cb25 |
c4 = GSL_MIN(1.0, (ascale / safemin) * bsize);
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
c4 = 1.0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (ascale <= 1.0 || bsize <= 1.0)
|
|
Packit |
67cb25 |
c5 = GSL_MIN(1.0, ascale * bsize);
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
c5 = 1.0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* scale first eigenvalue */
|
|
Packit |
67cb25 |
wabs = fabs(*wr1) + fabs(*wi);
|
|
Packit |
67cb25 |
wsize = GSL_MAX(safemin,
|
|
Packit |
67cb25 |
GSL_MAX(c1,
|
|
Packit |
67cb25 |
GSL_MAX(fuzzy1 * (wabs*c2 + c3),
|
|
Packit |
67cb25 |
GSL_MIN(c4, 0.5 * GSL_MAX(wabs, c5)))));
|
|
Packit |
67cb25 |
if (wsize != 1.0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
wscale = 1.0 / wsize;
|
|
Packit |
67cb25 |
if (wsize > 1.0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
*scale1 = (GSL_MAX(ascale, bsize) * wscale) *
|
|
Packit |
67cb25 |
GSL_MIN(ascale, bsize);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
*scale1 = (GSL_MIN(ascale, bsize) * wscale) *
|
|
Packit |
67cb25 |
GSL_MAX(ascale, bsize);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
*wr1 *= wscale;
|
|
Packit |
67cb25 |
if (*wi != 0.0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
*wi *= wscale;
|
|
Packit |
67cb25 |
*wr2 = *wr1;
|
|
Packit |
67cb25 |
*scale2 = *scale1;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
*scale1 = ascale * bsize;
|
|
Packit |
67cb25 |
*scale2 = *scale1;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* scale second eigenvalue if real */
|
|
Packit |
67cb25 |
if (*wi == 0.0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
wsize = GSL_MAX(safemin,
|
|
Packit |
67cb25 |
GSL_MAX(c1,
|
|
Packit |
67cb25 |
GSL_MAX(fuzzy1 * (fabs(*wr2) * c2 + c3),
|
|
Packit |
67cb25 |
GSL_MIN(c4, 0.5 * GSL_MAX(fabs(*wr2), c5)))));
|
|
Packit |
67cb25 |
if (wsize != 1.0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
wscale = 1.0 / wsize;
|
|
Packit |
67cb25 |
if (wsize > 1.0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
*scale2 = (GSL_MAX(ascale, bsize) * wscale) *
|
|
Packit |
67cb25 |
GSL_MIN(ascale, bsize);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
*scale2 = (GSL_MIN(ascale, bsize) * wscale) *
|
|
Packit |
67cb25 |
GSL_MAX(ascale, bsize);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
*wr2 *= wscale;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
*scale2 = ascale * bsize;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
} /* gsl_schur_gen_eigvals() */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
gsl_schur_solve_equation()
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Solve the equation which comes up in the back substitution
|
|
Packit |
67cb25 |
when computing eigenvectors corresponding to real eigenvalues.
|
|
Packit |
67cb25 |
The equation that is solved is:
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
(ca*A - z*D)*x = s*b
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
where
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
A is n-by-n with n = 1 or 2
|
|
Packit |
67cb25 |
D is a n-by-n diagonal matrix
|
|
Packit |
67cb25 |
b and x are n-by-1 real vectors
|
|
Packit |
67cb25 |
s is a scaling factor set by this function to prevent overflow in x
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Inputs: ca - coefficient multiplying A
|
|
Packit |
67cb25 |
A - square matrix (n-by-n)
|
|
Packit |
67cb25 |
z - real scalar (eigenvalue)
|
|
Packit |
67cb25 |
d1 - (1,1) element in diagonal matrix D
|
|
Packit |
67cb25 |
d2 - (2,2) element in diagonal matrix D
|
|
Packit |
67cb25 |
b - right hand side vector
|
|
Packit |
67cb25 |
x - (output) where to store solution
|
|
Packit |
67cb25 |
s - (output) scale factor
|
|
Packit |
67cb25 |
xnorm - (output) infinity norm of X
|
|
Packit |
67cb25 |
smin - lower bound on singular values of A - if ca*A - z*D
|
|
Packit |
67cb25 |
is less than this value, we'll use smin*I instead.
|
|
Packit |
67cb25 |
This value should be a safe distance above underflow.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Return: success
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Notes: 1) A and b are not changed on output
|
|
Packit |
67cb25 |
2) Based on lapack routine DLALN2
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_schur_solve_equation(double ca, const gsl_matrix *A, double z,
|
|
Packit |
67cb25 |
double d1, double d2, const gsl_vector *b,
|
|
Packit |
67cb25 |
gsl_vector *x, double *s, double *xnorm,
|
|
Packit |
67cb25 |
double smin)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
size_t N = A->size1;
|
|
Packit |
67cb25 |
double bnorm;
|
|
Packit |
67cb25 |
double scale = 1.0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (N == 1)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double c, /* denominator */
|
|
Packit |
67cb25 |
cnorm; /* |c| */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* we have a 1-by-1 (real) scalar system to solve */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
c = ca * gsl_matrix_get(A, 0, 0) - z * d1;
|
|
Packit |
67cb25 |
cnorm = fabs(c);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (cnorm < smin)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* set c = smin*I */
|
|
Packit |
67cb25 |
c = smin;
|
|
Packit |
67cb25 |
cnorm = smin;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* check scaling for x = b / c */
|
|
Packit |
67cb25 |
bnorm = fabs(gsl_vector_get(b, 0));
|
|
Packit |
67cb25 |
if (cnorm < 1.0 && bnorm > 1.0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if (bnorm > GSL_SCHUR_BIGNUM*cnorm)
|
|
Packit |
67cb25 |
scale = 1.0 / bnorm;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* compute x */
|
|
Packit |
67cb25 |
gsl_vector_set(x, 0, gsl_vector_get(b, 0) * scale / c);
|
|
Packit |
67cb25 |
*xnorm = fabs(gsl_vector_get(x, 0));
|
|
Packit |
67cb25 |
} /* if (N == 1) */
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double cr[2][2];
|
|
Packit |
67cb25 |
double *crv;
|
|
Packit |
67cb25 |
double cmax;
|
|
Packit |
67cb25 |
size_t icmax, j;
|
|
Packit |
67cb25 |
double bval1, bval2;
|
|
Packit |
67cb25 |
double ur11, ur12, ur22, ur11r;
|
|
Packit |
67cb25 |
double cr21, cr22;
|
|
Packit |
67cb25 |
double lr21;
|
|
Packit |
67cb25 |
double b1, b2, bbnd;
|
|
Packit |
67cb25 |
double x1, x2;
|
|
Packit |
67cb25 |
double temp;
|
|
Packit |
67cb25 |
size_t ipivot[4][4] = { { 0, 1, 2, 3 },
|
|
Packit |
67cb25 |
{ 1, 0, 3, 2 },
|
|
Packit |
67cb25 |
{ 2, 3, 0, 1 },
|
|
Packit |
67cb25 |
{ 3, 2, 1, 0 } };
|
|
Packit |
67cb25 |
int rswap[4] = { 0, 1, 0, 1 };
|
|
Packit |
67cb25 |
int zswap[4] = { 0, 0, 1, 1 };
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
* we have a 2-by-2 real system to solve:
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* ( ca * [ A11 A12 ] - z * [ D1 0 ] ) [ x1 ] = [ b1 ]
|
|
Packit |
67cb25 |
* ( [ A21 A22 ] [ 0 D2 ] ) [ x2 ] [ b2 ]
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* (z real)
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
crv = (double *) cr;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
* compute the real part of C = ca*A - z*D - use column ordering
|
|
Packit |
67cb25 |
* here since porting from lapack
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
cr[0][0] = ca * gsl_matrix_get(A, 0, 0) - z * d1;
|
|
Packit |
67cb25 |
cr[1][1] = ca * gsl_matrix_get(A, 1, 1) - z * d2;
|
|
Packit |
67cb25 |
cr[0][1] = ca * gsl_matrix_get(A, 1, 0);
|
|
Packit |
67cb25 |
cr[1][0] = ca * gsl_matrix_get(A, 0, 1);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* find the largest element in C */
|
|
Packit |
67cb25 |
cmax = 0.0;
|
|
Packit |
67cb25 |
icmax = 0;
|
|
Packit |
67cb25 |
for (j = 0; j < 4; ++j)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if (fabs(crv[j]) > cmax)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
cmax = fabs(crv[j]);
|
|
Packit |
67cb25 |
icmax = j;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
bval1 = gsl_vector_get(b, 0);
|
|
Packit |
67cb25 |
bval2 = gsl_vector_get(b, 1);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* if norm(C) < smin, use smin*I */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (cmax < smin)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
bnorm = GSL_MAX(fabs(bval1), fabs(bval2));
|
|
Packit |
67cb25 |
if (smin < 1.0 && bnorm > 1.0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if (bnorm > GSL_SCHUR_BIGNUM*smin)
|
|
Packit |
67cb25 |
scale = 1.0 / bnorm;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
temp = scale / smin;
|
|
Packit |
67cb25 |
gsl_vector_set(x, 0, temp * bval1);
|
|
Packit |
67cb25 |
gsl_vector_set(x, 1, temp * bval2);
|
|
Packit |
67cb25 |
*xnorm = temp * bnorm;
|
|
Packit |
67cb25 |
*s = scale;
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* gaussian elimination with complete pivoting */
|
|
Packit |
67cb25 |
ur11 = crv[icmax];
|
|
Packit |
67cb25 |
cr21 = crv[ipivot[1][icmax]];
|
|
Packit |
67cb25 |
ur12 = crv[ipivot[2][icmax]];
|
|
Packit |
67cb25 |
cr22 = crv[ipivot[3][icmax]];
|
|
Packit |
67cb25 |
ur11r = 1.0 / ur11;
|
|
Packit |
67cb25 |
lr21 = ur11r * cr21;
|
|
Packit |
67cb25 |
ur22 = cr22 - ur12 * lr21;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* if smaller pivot < smin, use smin */
|
|
Packit |
67cb25 |
if (fabs(ur22) < smin)
|
|
Packit |
67cb25 |
ur22 = smin;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (rswap[icmax])
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
b1 = bval2;
|
|
Packit |
67cb25 |
b2 = bval1;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
b1 = bval1;
|
|
Packit |
67cb25 |
b2 = bval2;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
b2 -= lr21 * b1;
|
|
Packit |
67cb25 |
bbnd = GSL_MAX(fabs(b1 * (ur22 * ur11r)), fabs(b2));
|
|
Packit |
67cb25 |
if (bbnd > 1.0 && fabs(ur22) < 1.0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if (bbnd >= GSL_SCHUR_BIGNUM * fabs(ur22))
|
|
Packit |
67cb25 |
scale = 1.0 / bbnd;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
x2 = (b2 * scale) / ur22;
|
|
Packit |
67cb25 |
x1 = (scale * b1) * ur11r - x2 * (ur11r * ur12);
|
|
Packit |
67cb25 |
if (zswap[icmax])
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_vector_set(x, 0, x2);
|
|
Packit |
67cb25 |
gsl_vector_set(x, 1, x1);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_vector_set(x, 0, x1);
|
|
Packit |
67cb25 |
gsl_vector_set(x, 1, x2);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
*xnorm = GSL_MAX(fabs(x1), fabs(x2));
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* further scaling if norm(A) norm(X) > overflow */
|
|
Packit |
67cb25 |
if (*xnorm > 1.0 && cmax > 1.0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if (*xnorm > GSL_SCHUR_BIGNUM / cmax)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
temp = cmax / GSL_SCHUR_BIGNUM;
|
|
Packit |
67cb25 |
gsl_blas_dscal(temp, x);
|
|
Packit |
67cb25 |
*xnorm *= temp;
|
|
Packit |
67cb25 |
scale *= temp;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
} /* if (N == 2) */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
*s = scale;
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
} /* gsl_schur_solve_equation() */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
gsl_schur_solve_equation_z()
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Solve the equation which comes up in the back substitution
|
|
Packit |
67cb25 |
when computing eigenvectors corresponding to complex eigenvalues.
|
|
Packit |
67cb25 |
The equation that is solved is:
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
(ca*A - z*D)*x = s*b
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
where
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
A is n-by-n with n = 1 or 2
|
|
Packit |
67cb25 |
D is a n-by-n diagonal matrix
|
|
Packit |
67cb25 |
b and x are n-by-1 complex vectors
|
|
Packit |
67cb25 |
s is a scaling factor set by this function to prevent overflow in x
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Inputs: ca - coefficient multiplying A
|
|
Packit |
67cb25 |
A - square matrix (n-by-n)
|
|
Packit |
67cb25 |
z - complex scalar (eigenvalue)
|
|
Packit |
67cb25 |
d1 - (1,1) element in diagonal matrix D
|
|
Packit |
67cb25 |
d2 - (2,2) element in diagonal matrix D
|
|
Packit |
67cb25 |
b - right hand side vector
|
|
Packit |
67cb25 |
x - (output) where to store solution
|
|
Packit |
67cb25 |
s - (output) scale factor
|
|
Packit |
67cb25 |
xnorm - (output) infinity norm of X
|
|
Packit |
67cb25 |
smin - lower bound on singular values of A - if ca*A - z*D
|
|
Packit |
67cb25 |
is less than this value, we'll use smin*I instead.
|
|
Packit |
67cb25 |
This value should be a safe distance above underflow.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Notes: 1) A and b are not changed on output
|
|
Packit |
67cb25 |
2) Based on lapack routine DLALN2
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_schur_solve_equation_z(double ca, const gsl_matrix *A, gsl_complex *z,
|
|
Packit |
67cb25 |
double d1, double d2,
|
|
Packit |
67cb25 |
const gsl_vector_complex *b,
|
|
Packit |
67cb25 |
gsl_vector_complex *x, double *s, double *xnorm,
|
|
Packit |
67cb25 |
double smin)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
size_t N = A->size1;
|
|
Packit |
67cb25 |
double scale = 1.0;
|
|
Packit |
67cb25 |
double bnorm;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (N == 1)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double cr, /* denominator */
|
|
Packit |
67cb25 |
ci,
|
|
Packit |
67cb25 |
cnorm; /* |c| */
|
|
Packit |
67cb25 |
gsl_complex bval, c, xval, tmp;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* we have a 1-by-1 (complex) scalar system to solve */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* c = ca*a - z*d1 */
|
|
Packit |
67cb25 |
cr = ca * gsl_matrix_get(A, 0, 0) - GSL_REAL(*z) * d1;
|
|
Packit |
67cb25 |
ci = -GSL_IMAG(*z) * d1;
|
|
Packit |
67cb25 |
cnorm = fabs(cr) + fabs(ci);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (cnorm < smin)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* set c = smin*I */
|
|
Packit |
67cb25 |
cr = smin;
|
|
Packit |
67cb25 |
ci = 0.0;
|
|
Packit |
67cb25 |
cnorm = smin;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* check scaling for x = b / c */
|
|
Packit |
67cb25 |
bval = gsl_vector_complex_get(b, 0);
|
|
Packit |
67cb25 |
bnorm = fabs(GSL_REAL(bval)) + fabs(GSL_IMAG(bval));
|
|
Packit |
67cb25 |
if (cnorm < 1.0 && bnorm > 1.0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if (bnorm > GSL_SCHUR_BIGNUM*cnorm)
|
|
Packit |
67cb25 |
scale = 1.0 / bnorm;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* compute x */
|
|
Packit |
67cb25 |
GSL_SET_COMPLEX(&tmp, scale*GSL_REAL(bval), scale*GSL_IMAG(bval));
|
|
Packit |
67cb25 |
GSL_SET_COMPLEX(&c, cr, ci);
|
|
Packit |
67cb25 |
xval = gsl_complex_div(tmp, c);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_vector_complex_set(x, 0, xval);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
*xnorm = fabs(GSL_REAL(xval)) + fabs(GSL_IMAG(xval));
|
|
Packit |
67cb25 |
} /* if (N == 1) */
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double cr[2][2], ci[2][2];
|
|
Packit |
67cb25 |
double *civ, *crv;
|
|
Packit |
67cb25 |
double cmax;
|
|
Packit |
67cb25 |
gsl_complex bval1, bval2;
|
|
Packit |
67cb25 |
gsl_complex xval1, xval2;
|
|
Packit |
67cb25 |
double xr1, xi1;
|
|
Packit |
67cb25 |
size_t icmax;
|
|
Packit |
67cb25 |
size_t j;
|
|
Packit |
67cb25 |
double temp;
|
|
Packit |
67cb25 |
double ur11, ur12, ur22, ui11, ui12, ui22, ur11r, ui11r;
|
|
Packit |
67cb25 |
double ur12s, ui12s;
|
|
Packit |
67cb25 |
double u22abs;
|
|
Packit |
67cb25 |
double lr21, li21;
|
|
Packit |
67cb25 |
double cr21, cr22, ci21, ci22;
|
|
Packit |
67cb25 |
double br1, bi1, br2, bi2, bbnd;
|
|
Packit |
67cb25 |
gsl_complex b1, b2;
|
|
Packit |
67cb25 |
size_t ipivot[4][4] = { { 0, 1, 2, 3 },
|
|
Packit |
67cb25 |
{ 1, 0, 3, 2 },
|
|
Packit |
67cb25 |
{ 2, 3, 0, 1 },
|
|
Packit |
67cb25 |
{ 3, 2, 1, 0 } };
|
|
Packit |
67cb25 |
int rswap[4] = { 0, 1, 0, 1 };
|
|
Packit |
67cb25 |
int zswap[4] = { 0, 0, 1, 1 };
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
* complex 2-by-2 system:
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* ( ca * [ A11 A12 ] - z * [ D1 0 ] ) [ X1 ] = [ B1 ]
|
|
Packit |
67cb25 |
* ( [ A21 A22 ] [ 0 D2] ) [ X2 ] [ B2 ]
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* (z complex)
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* where the X and B values are complex.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
civ = (double *) ci;
|
|
Packit |
67cb25 |
crv = (double *) cr;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
* compute the real part of C = ca*A - z*D - use column ordering
|
|
Packit |
67cb25 |
* here since porting from lapack
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
cr[0][0] = ca*gsl_matrix_get(A, 0, 0) - GSL_REAL(*z)*d1;
|
|
Packit |
67cb25 |
cr[1][1] = ca*gsl_matrix_get(A, 1, 1) - GSL_REAL(*z)*d2;
|
|
Packit |
67cb25 |
cr[0][1] = ca*gsl_matrix_get(A, 1, 0);
|
|
Packit |
67cb25 |
cr[1][0] = ca*gsl_matrix_get(A, 0, 1);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* compute the imaginary part */
|
|
Packit |
67cb25 |
ci[0][0] = -GSL_IMAG(*z) * d1;
|
|
Packit |
67cb25 |
ci[0][1] = 0.0;
|
|
Packit |
67cb25 |
ci[1][0] = 0.0;
|
|
Packit |
67cb25 |
ci[1][1] = -GSL_IMAG(*z) * d2;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
cmax = 0.0;
|
|
Packit |
67cb25 |
icmax = 0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (j = 0; j < 4; ++j)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if (fabs(crv[j]) + fabs(civ[j]) > cmax)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
cmax = fabs(crv[j]) + fabs(civ[j]);
|
|
Packit |
67cb25 |
icmax = j;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
bval1 = gsl_vector_complex_get(b, 0);
|
|
Packit |
67cb25 |
bval2 = gsl_vector_complex_get(b, 1);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* if norm(C) < smin, use smin*I */
|
|
Packit |
67cb25 |
if (cmax < smin)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
bnorm = GSL_MAX(fabs(GSL_REAL(bval1)) + fabs(GSL_IMAG(bval1)),
|
|
Packit |
67cb25 |
fabs(GSL_REAL(bval2)) + fabs(GSL_IMAG(bval2)));
|
|
Packit |
67cb25 |
if (smin < 1.0 && bnorm > 1.0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if (bnorm > GSL_SCHUR_BIGNUM*smin)
|
|
Packit |
67cb25 |
scale = 1.0 / bnorm;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
temp = scale / smin;
|
|
Packit |
67cb25 |
xval1 = gsl_complex_mul_real(bval1, temp);
|
|
Packit |
67cb25 |
xval2 = gsl_complex_mul_real(bval2, temp);
|
|
Packit |
67cb25 |
gsl_vector_complex_set(x, 0, xval1);
|
|
Packit |
67cb25 |
gsl_vector_complex_set(x, 1, xval2);
|
|
Packit |
67cb25 |
*xnorm = temp * bnorm;
|
|
Packit |
67cb25 |
*s = scale;
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* gaussian elimination with complete pivoting */
|
|
Packit |
67cb25 |
ur11 = crv[icmax];
|
|
Packit |
67cb25 |
ui11 = civ[icmax];
|
|
Packit |
67cb25 |
cr21 = crv[ipivot[1][icmax]];
|
|
Packit |
67cb25 |
ci21 = civ[ipivot[1][icmax]];
|
|
Packit |
67cb25 |
ur12 = crv[ipivot[2][icmax]];
|
|
Packit |
67cb25 |
ui12 = civ[ipivot[2][icmax]];
|
|
Packit |
67cb25 |
cr22 = crv[ipivot[3][icmax]];
|
|
Packit |
67cb25 |
ci22 = civ[ipivot[3][icmax]];
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (icmax == 0 || icmax == 3)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* off diagonals of pivoted C are real */
|
|
Packit |
67cb25 |
if (fabs(ur11) > fabs(ui11))
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
temp = ui11 / ur11;
|
|
Packit |
67cb25 |
ur11r = 1.0 / (ur11 * (1.0 + temp*temp));
|
|
Packit |
67cb25 |
ui11r = -temp * ur11r;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
temp = ur11 / ui11;
|
|
Packit |
67cb25 |
ui11r = -1.0 / (ui11 * (1.0 + temp*temp));
|
|
Packit |
67cb25 |
ur11r = -temp*ui11r;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
lr21 = cr21 * ur11r;
|
|
Packit |
67cb25 |
li21 = cr21 * ui11r;
|
|
Packit |
67cb25 |
ur12s = ur12 * ur11r;
|
|
Packit |
67cb25 |
ui12s = ur12 * ui11r;
|
|
Packit |
67cb25 |
ur22 = cr22 - ur12 * lr21;
|
|
Packit |
67cb25 |
ui22 = ci22 - ur12 * li21;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* diagonals of pivoted C are real */
|
|
Packit |
67cb25 |
ur11r = 1.0 / ur11;
|
|
Packit |
67cb25 |
ui11r = 0.0;
|
|
Packit |
67cb25 |
lr21 = cr21 * ur11r;
|
|
Packit |
67cb25 |
li21 = ci21 * ur11r;
|
|
Packit |
67cb25 |
ur12s = ur12 * ur11r;
|
|
Packit |
67cb25 |
ui12s = ui12 * ur11r;
|
|
Packit |
67cb25 |
ur22 = cr22 - ur12 * lr21 + ui12 * li21;
|
|
Packit |
67cb25 |
ui22 = -ur12 * li21 - ui12 * lr21;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
u22abs = fabs(ur22) + fabs(ui22);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* if smaller pivot < smin, use smin */
|
|
Packit |
67cb25 |
if (u22abs < smin)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
ur22 = smin;
|
|
Packit |
67cb25 |
ui22 = 0.0;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (rswap[icmax])
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
br2 = GSL_REAL(bval1);
|
|
Packit |
67cb25 |
bi2 = GSL_IMAG(bval1);
|
|
Packit |
67cb25 |
br1 = GSL_REAL(bval2);
|
|
Packit |
67cb25 |
bi1 = GSL_IMAG(bval2);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
br1 = GSL_REAL(bval1);
|
|
Packit |
67cb25 |
bi1 = GSL_IMAG(bval1);
|
|
Packit |
67cb25 |
br2 = GSL_REAL(bval2);
|
|
Packit |
67cb25 |
bi2 = GSL_IMAG(bval2);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
br2 += li21*bi1 - lr21*br1;
|
|
Packit |
67cb25 |
bi2 -= li21*br1 + lr21*bi1;
|
|
Packit |
67cb25 |
bbnd = GSL_MAX((fabs(br1) + fabs(bi1)) *
|
|
Packit |
67cb25 |
(u22abs * (fabs(ur11r) + fabs(ui11r))),
|
|
Packit |
67cb25 |
fabs(br2) + fabs(bi2));
|
|
Packit |
67cb25 |
if (bbnd > 1.0 && u22abs < 1.0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if (bbnd >= GSL_SCHUR_BIGNUM*u22abs)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
scale = 1.0 / bbnd;
|
|
Packit |
67cb25 |
br1 *= scale;
|
|
Packit |
67cb25 |
bi1 *= scale;
|
|
Packit |
67cb25 |
br2 *= scale;
|
|
Packit |
67cb25 |
bi2 *= scale;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
GSL_SET_COMPLEX(&b1, br2, bi2);
|
|
Packit |
67cb25 |
GSL_SET_COMPLEX(&b2, ur22, ui22);
|
|
Packit |
67cb25 |
xval2 = gsl_complex_div(b1, b2);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
xr1 = ur11r*br1 - ui11r*bi1 - ur12s*GSL_REAL(xval2) + ui12s*GSL_IMAG(xval2);
|
|
Packit |
67cb25 |
xi1 = ui11r*br1 + ur11r*bi1 - ui12s*GSL_REAL(xval2) - ur12s*GSL_IMAG(xval2);
|
|
Packit |
67cb25 |
GSL_SET_COMPLEX(&xval1, xr1, xi1);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (zswap[icmax])
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_vector_complex_set(x, 0, xval2);
|
|
Packit |
67cb25 |
gsl_vector_complex_set(x, 1, xval1);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_vector_complex_set(x, 0, xval1);
|
|
Packit |
67cb25 |
gsl_vector_complex_set(x, 1, xval2);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
*xnorm = GSL_MAX(fabs(GSL_REAL(xval1)) + fabs(GSL_IMAG(xval1)),
|
|
Packit |
67cb25 |
fabs(GSL_REAL(xval2)) + fabs(GSL_IMAG(xval2)));
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* further scaling if norm(A) norm(X) > overflow */
|
|
Packit |
67cb25 |
if (*xnorm > 1.0 && cmax > 1.0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if (*xnorm > GSL_SCHUR_BIGNUM / cmax)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
temp = cmax / GSL_SCHUR_BIGNUM;
|
|
Packit |
67cb25 |
gsl_blas_zdscal(temp, x);
|
|
Packit |
67cb25 |
*xnorm *= temp;
|
|
Packit |
67cb25 |
scale *= temp;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
} /* if (N == 2) */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
*s = scale;
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
} /* gsl_schur_solve_equation_z() */
|