cern.colt.matrix.tfloat.algo.solver.FloatHyBR Maven / Gradle / Ivy
Go to download
Show more of this group Show more artifacts with this name
Show all versions of parallelcolt Show documentation
Show all versions of parallelcolt Show documentation
Parallel Colt is a multithreaded version of Colt - a library for high performance scientific computing in Java. It contains efficient algorithms for data analysis, linear algebra, multi-dimensional arrays, Fourier transforms, statistics and histogramming.
The newest version!
/*
* Copyright (C) 2009 Piotr Wendykier
*
* This program is free software; you can redistribute it and/or modify it
* under the terms of the GNU Lesser General Public License as published by the
* Free Software Foundation; either version 2.1 of the License, or (at your
* option) any later version.
*
* This library 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 Lesser General Public License
* for more details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with this library; if not, write to the Free Software Foundation,
* Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*/
package cern.colt.matrix.tfloat.algo.solver;
import optimization.FloatFmin;
import optimization.FloatFmin_methods;
import cern.colt.list.tfloat.FloatArrayList;
import cern.colt.matrix.tdouble.algo.solver.HyBRInnerSolver;
import cern.colt.matrix.tdouble.algo.solver.HyBRRegularizationMethod;
import cern.colt.matrix.tfloat.FloatFactory2D;
import cern.colt.matrix.tfloat.FloatMatrix1D;
import cern.colt.matrix.tfloat.FloatMatrix2D;
import cern.colt.matrix.tfloat.algo.DenseFloatAlgebra;
import cern.colt.matrix.tfloat.algo.decomposition.DenseFloatSingularValueDecomposition;
import cern.colt.matrix.tfloat.algo.solver.preconditioner.FloatIdentity;
import cern.colt.matrix.tfloat.algo.solver.preconditioner.FloatPreconditioner;
import cern.colt.matrix.tfloat.impl.DenseColumnFloatMatrix2D;
import cern.colt.matrix.tfloat.impl.DenseFloatMatrix1D;
import cern.colt.matrix.tfloat.impl.DenseFloatMatrix2D;
import cern.jet.math.tfloat.FloatFunctions;
import cern.jet.stat.tfloat.FloatDescriptive;
/**
* HyBR is a Hybrid Bidiagonalization Regularization method used for solving
* large-scale, ill-posed inverse problems of the form: b = A*x + noise The
* method combines an iterative Lanczos Bidiagonalization (LBD) Method with an
* SVD-based regularization method to stabilize the semiconvergence behavior
* that is characteristic of many ill-posed problems. The code is derived from
* RestoreTools: An Object Oriented Matlab Package for Image Restoration written
* by James G. Nagy and several of his students, including Julianne Chung,
* Katrina Palmer, Lisa Perrone, and Ryan Wright.
*
*
* References:
*
* [1] Paige and Saunders, "LSQR an algorithm for sparse linear equations an
* sparse least squares", ACM Trans. Math Software, 8 (1982), pp. 43-71.
*
*
* [2] Bjorck, Grimme and Van Dooren, "An implicit shift bidiagonalization
* algorithm for ill-posed systems", BIT 34 (11994), pp. 520-534.
*
*
* [3] Chung, Nagy and O'Leary, "A Weighted GCV Method for Lanczos Hybrid
* Regularization", Elec. Trans. Numer. Anal., 28 (2008), pp. 149--167.
*
*
*
* @author Piotr Wendykier ([email protected])
*
*/
public class FloatHyBR extends AbstractFloatIterativeSolver {
private HyBRInnerSolver innerSolver;
private HyBRRegularizationMethod regMethod;
private float regPar;
private float omega;
private boolean reorth;
private int begReg;
private float flatTol;
private boolean computeRnrm;
private static final DenseFloatAlgebra alg = DenseFloatAlgebra.DEFAULT;
private static final float FMIN_TOL = 1.0e-4f;
/**
* Creates new instance of HyBR solver with default parameters:
*
* innerSolver = HyBR.InnerSolver.TIKHONOV
* regularizationMethod = HyBR.RegularizationMethod.ADAPTWGCV
* regularizationParameter = 0
* omega = 0
* reorthogonalize = false
* beginRegularization = 2
* flatTolerance = 1e-6
* computeRnrm = false;
*/
public FloatHyBR() {
this(HyBRInnerSolver.TIKHONOV, HyBRRegularizationMethod.ADAPTWGCV, 0, 0, false, 2, 1e-4f, false);
}
/**
* Creates new instance of HyBR solver.
*
* @param innerSolver
* solver for the inner problem
* @param regularizationMethod
* a method for choosing a regularization parameter
* @param regularizationParameter
* if regularizationMethod == HyBR.RegularizationMethod.NONE then
* the regularization parameter has to be specified here (value
* from the interval (0,1))
* @param omega
* regularizationMethod == HyBR.RegularizationMethod.WGCV then
* omega has to be specified here (must be nonnegative)
* @param reorthogonalize
* if thue then Lanczos subspaces are reorthogonalized
* @param beginRegularization
* begin regularization after this iteration (must be at least 2)
* @param flatTolerance
* tolerance for detecting flatness in the GCV curve as a
* stopping criteria (must be nonnegative)
*@param computeRnrm
* if true then the norm of relative residual is computed
*/
public FloatHyBR(HyBRInnerSolver innerSolver, HyBRRegularizationMethod regularizationMethod,
float regularizationParameter, float omega, boolean reorthogonalize, int beginRegularization,
float flatTolerance, boolean computeRnrm) {
this.innerSolver = innerSolver;
this.regMethod = regularizationMethod;
if ((regularizationParameter < 0.0) || (regularizationParameter > 1.0)) {
throw new IllegalArgumentException("regularizationParameter must be a number between 0 and 1.");
}
this.regPar = regularizationParameter;
if (omega < 0.0) {
throw new IllegalArgumentException("omega must be a nonnegative number.");
}
this.omega = omega;
this.reorth = reorthogonalize;
if (beginRegularization < 2) {
throw new IllegalArgumentException("beginRegularization must be greater or equal 2");
}
this.begReg = beginRegularization;
if (flatTolerance < 0.0) {
throw new IllegalArgumentException("flatTolerance must be a nonnegative number.");
}
this.flatTol = flatTolerance;
this.computeRnrm = computeRnrm;
this.iter = new HyBRFloatIterationMonitor();
}
public FloatMatrix1D solve(FloatMatrix2D A, FloatMatrix1D b, FloatMatrix1D x)
throws IterativeSolverFloatNotConvergedException {
if (!(iter instanceof HyBRFloatIterationMonitor)) {
this.iter = new HyBRFloatIterationMonitor();
}
checkSizes(A, b, x);
int rows = A.rows();
int columns = A.columns();
boolean bump = false;
boolean warning = false;
float rnrm = -1.0f;
int iterationsSave = 0;
float alpha = 0;
float alphaSave = 0;
float beta = 0;
HyBRInnerSolver inSolver = HyBRInnerSolver.NONE;
FloatLBD lbd;
FloatMatrix1D v;
FloatMatrix1D work;
FloatMatrix2D Ub, Vb;
FloatMatrix1D f = null;
FloatMatrix1D xSave = null;
float[] sv;
FloatArrayList omegaList = new FloatArrayList(new float[begReg - 2]);
FloatArrayList GCV = new FloatArrayList(new float[begReg - 2]);
FloatMatrix2D U = new DenseFloatMatrix2D(1, (int) b.size());
FloatMatrix2D C = null;
FloatMatrix2D V = null;
DenseFloatSingularValueDecomposition svd;
if (computeRnrm) {
work = b.copy();
A.zMult(x, work, -1, 1, false);
rnrm = alg.norm2(work);
}
if (M instanceof FloatIdentity) {
beta = alg.norm2(b);
U.viewRow(0).assign(b, FloatFunctions.multSecond(1.0f / beta));
lbd = new FloatSimpleLBD(A, U, reorth);
} else {
work = new DenseFloatMatrix1D((int) b.size());
work = M.apply(b, work);
beta = alg.norm2(work);
U.viewRow(0).assign(work, FloatFunctions.multSecond(1.0f / beta));
lbd = new FloatPLBD(M, A, U, reorth);
}
for (iter.setFirst(); !iter.converged(rnrm, x); iter.next()) {
lbd.apply();
U = lbd.getU();
C = lbd.getC();
V = lbd.getV();
v = new DenseFloatMatrix1D(C.columns() + 1);
v.setQuick(0, beta);
int i = iter.iterations();
if (i >= 1) {
if (i >= begReg - 1) {
inSolver = innerSolver;
}
switch (inSolver) {
case TIKHONOV:
svd = alg.svd(C);
Ub = svd.getU();
sv = svd.getSingularValues();
Vb = svd.getV();
if (regMethod == HyBRRegularizationMethod.ADAPTWGCV) {
work = new DenseFloatMatrix1D(Ub.rows());
Ub.zMult(v, work, 1, 0, true);
omegaList.add(Math.min(1, findOmega(work, sv)));
omega = FloatDescriptive.mean(omegaList);
}
f = new DenseFloatMatrix1D(Vb.rows());
alpha = tikhonovSolver(Ub, sv, Vb, v, f);
((HyBRFloatIterationMonitor) iter).setRegularizationParameter(alpha);
GCV.add(GCVstopfun(alpha, Ub.viewRow(0), sv, beta, rows, columns));
if (i > 1) {
if (Math.abs((GCV.getQuick(i - 1) - GCV.getQuick(i - 2))) / GCV.get(begReg - 2) < flatTol) {
V.zMult(f, x);
((HyBRFloatIterationMonitor) iter)
.setStoppingCondition(HyBRFloatIterationMonitor.HyBRStoppingCondition.FLAT_GCV_CURVE);
if (computeRnrm) {
work = b.copy();
A.zMult(x, work, -1, 1, false);
((HyBRFloatIterationMonitor) iter).residual = alg.norm2(work);
}
return x;
} else if ((warning == true) && (GCV.size() > iterationsSave + 3)) {
for (int j = iterationsSave; j < GCV.size(); j++) {
if (GCV.getQuick(iterationsSave - 1) > GCV.get(j)) {
bump = true;
}
}
if (bump == false) {
x.assign(xSave);
((HyBRFloatIterationMonitor) iter)
.setStoppingCondition(HyBRFloatIterationMonitor.HyBRStoppingCondition.MIN_OF_GCV_CURVE_WITHIN_WINDOW_OF_4_ITERATIONS);
((HyBRFloatIterationMonitor) iter).iter = iterationsSave;
if (computeRnrm) {
work = b.copy();
A.zMult(x, work, -1, 1, false);
((HyBRFloatIterationMonitor) iter).residual = alg.norm2(work);
}
((HyBRFloatIterationMonitor) iter).setRegularizationParameter(alphaSave);
return x;
} else {
bump = false;
warning = false;
iterationsSave = iter.getMaxIterations();
}
} else if (warning == false) {
if (GCV.get(i - 2) < GCV.get(i - 1)) {
warning = true;
xSave = new DenseFloatMatrix1D(V.rows());
alphaSave = alpha;
V.zMult(f, xSave);
iterationsSave = i;
}
}
}
break;
case NONE:
f = alg.solve(C, v);
break;
}
V.zMult(f, x);
if (computeRnrm) {
work = b.copy();
A.zMult(x, work, -1, 1, false);
rnrm = alg.norm2(work);
}
}
}
return x;
}
protected void checkSizes(FloatMatrix2D A, FloatMatrix1D b, FloatMatrix1D x) {
if (b.size() != A.rows())
throw new IllegalArgumentException("b.size() != A.rows()");
if (x.size() != A.columns())
throw new IllegalArgumentException("x.size() != A.columns()");
}
private float findOmega(FloatMatrix1D bhat, float[] s) {
int m = (int) bhat.size();
int n = s.length;
float alpha = s[n - 1];
float t0 = bhat.viewPart(n, m - n).aggregate(FloatFunctions.plus, FloatFunctions.square);
FloatMatrix1D s2 = new DenseFloatMatrix1D(s);
s2.assign(FloatFunctions.square);
float alpha2 = alpha * alpha;
FloatMatrix1D tt = s2.copy();
tt.assign(FloatFunctions.plus(alpha2));
tt.assign(FloatFunctions.inv);
float t1 = s2.aggregate(tt, FloatFunctions.plus, FloatFunctions.mult);
s2 = new DenseFloatMatrix1D(s);
s2.assign(FloatFunctions.mult(alpha));
s2.assign(bhat.viewPart(0, n), FloatFunctions.mult);
s2.assign(FloatFunctions.square);
FloatMatrix1D work = tt.copy();
work.assign(FloatFunctions.pow(3));
work.assign(FloatFunctions.abs);
float t3 = work.aggregate(s2, FloatFunctions.plus, FloatFunctions.mult);
work = new DenseFloatMatrix1D(s);
work.assign(tt, FloatFunctions.mult);
float t4 = work.aggregate(FloatFunctions.plus, FloatFunctions.square);
work = tt.copy();
work.assign(bhat.viewPart(0, n), FloatFunctions.mult);
work.assign(FloatFunctions.mult(alpha2));
float t5 = work.aggregate(FloatFunctions.plus, FloatFunctions.square);
s2 = new DenseFloatMatrix1D(s);
s2.assign(bhat.viewPart(0, n), FloatFunctions.mult);
s2.assign(FloatFunctions.square);
tt.assign(FloatFunctions.pow(3));
tt.assign(FloatFunctions.abs);
float v2 = tt.aggregate(s2, FloatFunctions.plus, FloatFunctions.mult);
return (m * alpha2 * v2) / (t1 * t3 + t4 * (t5 + t0));
}
private float tikhonovSolver(FloatMatrix2D U, float[] s, FloatMatrix2D V, FloatMatrix1D b, FloatMatrix1D x) {
TikFmin2D fmin;
FloatMatrix1D bhat = new DenseFloatMatrix1D(U.rows());
U.zMult(b, bhat, 1, 0, true);
float alpha = 0;
switch (regMethod) {
case GCV:
fmin = new TikFmin2D(bhat, s, 1);
alpha = FloatFmin.fmin(0, s[0], fmin, FMIN_TOL);
break;
case WGCV:
fmin = new TikFmin2D(bhat, s, omega);
alpha = FloatFmin.fmin(0, s[0], fmin, FMIN_TOL);
break;
case ADAPTWGCV:
fmin = new TikFmin2D(bhat, s, omega);
alpha = FloatFmin.fmin(0, s[0], fmin, FMIN_TOL);
break;
case NONE: // regularization parameter is given
alpha = regPar;
break;
}
FloatMatrix1D d = new DenseFloatMatrix1D(s);
d.assign(FloatFunctions.square);
d.assign(FloatFunctions.plus(alpha * alpha));
bhat = bhat.viewPart(0, s.length);
FloatMatrix1D S = new DenseFloatMatrix1D(s);
bhat.assign(S, FloatFunctions.mult);
bhat.assign(d, FloatFunctions.div);
V.zMult(bhat, x);
return alpha;
}
private static class TikFmin2D implements FloatFmin_methods {
FloatMatrix1D bhat;
float[] s;
float omega;
public TikFmin2D(FloatMatrix1D bhat, float[] s, float omega) {
this.bhat = bhat;
this.s = s;
this.omega = omega;
}
public float f_to_minimize(float alpha) {
int m = (int) bhat.size();
int n = s.length;
float t0 = bhat.viewPart(n, m - n).aggregate(FloatFunctions.plus, FloatFunctions.square);
FloatMatrix1D s2 = new DenseFloatMatrix1D(s);
s2.assign(FloatFunctions.square);
float alpha2 = alpha * alpha;
FloatMatrix1D work = s2.copy();
work.assign(FloatFunctions.plus(alpha2));
work.assign(FloatFunctions.inv);
FloatMatrix1D t1 = work.copy();
t1.assign(FloatFunctions.mult(alpha2));
FloatMatrix1D t2 = t1.copy();
t2.assign(bhat.viewPart(0, n), FloatFunctions.mult);
FloatMatrix1D t3 = work.copy();
t3.assign(s2, FloatFunctions.mult);
t3.assign(FloatFunctions.mult(1 - omega));
float denom = t3.aggregate(t1, FloatFunctions.plus, FloatFunctions.plus) + m - n;
return n * (t2.aggregate(FloatFunctions.plus, FloatFunctions.square) + t0) / (denom * denom);
}
}
private float GCVstopfun(float alpha, FloatMatrix1D u, float[] s, float beta, int rows, int columns) {
int k = s.length;
float beta2 = beta * beta;
FloatMatrix1D s2 = new DenseFloatMatrix1D(s);
s2.assign(FloatFunctions.square);
float alpha2 = alpha * alpha;
FloatMatrix1D t1 = s2.copy();
t1.assign(FloatFunctions.plus(alpha2));
t1.assign(FloatFunctions.inv);
FloatMatrix1D t2 = t1.copy();
t2.assign(u.viewPart(0, k), FloatFunctions.mult);
t2.assign(FloatFunctions.mult(alpha2));
float num = (float) (beta2
* (t2.aggregate(FloatFunctions.plus, FloatFunctions.square) + Math.pow(Math.abs(u.getQuick(k)), 2)) / columns);
float den = (rows - t1.aggregate(s2, FloatFunctions.plus, FloatFunctions.mult)) / columns;
den = den * den;
return num / den;
}
private interface FloatLBD {
public void apply();
public FloatMatrix2D getC();
public FloatMatrix2D getU();
public FloatMatrix2D getV();
}
private class FloatSimpleLBD implements FloatLBD {
private final DenseFloatAlgebra alg = DenseFloatAlgebra.DEFAULT;
private final FloatFactory2D factory = FloatFactory2D.dense;
private final FloatMatrix2D alphaBeta = new DenseFloatMatrix2D(2, 1);
private final FloatMatrix2D A;
private FloatMatrix2D C;
private FloatMatrix2D U;
private FloatMatrix2D V;
private boolean reorth;
private int counter = 1;
public FloatSimpleLBD(FloatMatrix2D A, FloatMatrix2D U, boolean reorth) {
this.A = A;
this.reorth = reorth;
this.U = U;
this.V = null;
this.C = null;
}
public void apply() {
if (reorth) {
int k = U.rows();
FloatMatrix1D u = null;
FloatMatrix1D v = null;
FloatMatrix1D column = null;
if (k == 1) {
v = A.zMult(U.viewRow(k - 1), v, 1, 0, true);
} else {
v = A.zMult(U.viewRow(k - 1), v, 1, 0, true);
column = V.viewColumn(k - 2);
v.assign(column, FloatFunctions.plusMultSecond(-C.getQuick(k - 1, k - 2)));
for (int j = 0; j < k - 1; j++) {
column = V.viewColumn(j);
v.assign(column, FloatFunctions.plusMultSecond(-column.zDotProduct(v)));
}
}
float alpha = alg.norm2(v);
v.assign(FloatFunctions.div(alpha));
u = A.zMult(v, u);
column = U.viewRow(k - 1);
u.assign(column, FloatFunctions.plusMultSecond(-alpha));
for (int j = 0; j < k; j++) {
column = U.viewRow(j);
u.assign(column, FloatFunctions.plusMultSecond(-column.zDotProduct(u)));
}
float beta = alg.norm2(u);
alphaBeta.setQuick(0, 0, alpha);
alphaBeta.setQuick(1, 0, beta);
u.assign(FloatFunctions.div(beta));
U = factory.appendRow(U, u);
if (V == null) {
V = new DenseColumnFloatMatrix2D((int) v.size(), 1);
V.assign((float[]) v.elements());
} else {
V = factory.appendColumn(V, v);
}
if (C == null) {
C = new DenseFloatMatrix2D(2, 1);
C.assign(alphaBeta);
} else {
C = factory.composeBidiagonal(C, alphaBeta);
}
} else {
FloatMatrix1D u = null;
FloatMatrix1D v = null;
FloatMatrix1D column = null;
if (counter == 1) {
v = A.zMult(U.viewRow(0), v, 1, 0, true);
} else {
v = A.zMult(U.viewRow(0), v, 1, 0, true);
column = V.viewColumn(counter - 2);
v.assign(column, FloatFunctions.plusMultSecond(-C.getQuick(counter - 1, counter - 2)));
}
float alpha = alg.norm2(v);
v.assign(FloatFunctions.div(alpha));
u = A.zMult(v, u);
column = U.viewRow(0);
u.assign(column, FloatFunctions.plusMultSecond(-alpha));
float beta = alg.norm2(u);
alphaBeta.setQuick(0, 0, alpha);
alphaBeta.setQuick(1, 0, beta);
u.assign(FloatFunctions.div(beta));
U.viewRow(0).assign(u);
if (V == null) {
V = new DenseColumnFloatMatrix2D((int) v.size(), 1);
V.assign((float[]) v.elements());
} else {
V = factory.appendColumn(V, v);
}
if (C == null) {
C = new DenseFloatMatrix2D(2, 1);
C.assign(alphaBeta);
} else {
C = factory.composeBidiagonal(C, alphaBeta);
}
counter++;
}
}
public FloatMatrix2D getC() {
return C;
}
public FloatMatrix2D getU() {
return U;
}
public FloatMatrix2D getV() {
return V;
}
}
private class FloatPLBD implements FloatLBD {
private final DenseFloatAlgebra alg = DenseFloatAlgebra.DEFAULT;
private final FloatFactory2D factory = FloatFactory2D.dense;
private final FloatMatrix2D alphaBeta = new DenseFloatMatrix2D(2, 1);
private final FloatPreconditioner M;
private final FloatMatrix2D A;
private FloatMatrix2D C;
private FloatMatrix2D U;
private FloatMatrix2D V;
private boolean reorth;
private int counter = 1;
public FloatPLBD(FloatPreconditioner M, FloatMatrix2D A, FloatMatrix2D U, boolean reorth) {
this.M = M;
this.A = A;
this.reorth = reorth;
this.U = U;
this.V = null;
this.C = null;
}
public void apply() {
if (reorth) {
int k = U.rows();
FloatMatrix1D u = null;
FloatMatrix1D v = null;
FloatMatrix1D row = null;
if (k == 1) {
row = U.viewRow(k - 1).copy();
row = M.transApply(row, row);
v = A.zMult(row, v, 1, 0, true);
} else {
row = U.viewRow(k - 1).copy();
row = M.transApply(row, row);
v = A.zMult(row, v, 1, 0, true);
row = V.viewColumn(k - 2);
v.assign(row, FloatFunctions.plusMultSecond(-C.getQuick(k - 1, k - 2)));
for (int j = 0; j < k - 1; j++) {
row = V.viewColumn(j);
v.assign(row, FloatFunctions.plusMultSecond(-row.zDotProduct(v)));
}
}
float alpha = alg.norm2(v);
v.assign(FloatFunctions.div(alpha));
row = A.zMult(v, row);
u = M.apply(row, u);
row = U.viewRow(k - 1);
u.assign(row, FloatFunctions.plusMultSecond(-alpha));
for (int j = 0; j < k; j++) {
row = U.viewRow(j);
u.assign(row, FloatFunctions.plusMultSecond(-row.zDotProduct(u)));
}
float beta = alg.norm2(u);
alphaBeta.setQuick(0, 0, alpha);
alphaBeta.setQuick(1, 0, beta);
u.assign(FloatFunctions.div(beta));
U = factory.appendRow(U, u);
if (V == null) {
V = new DenseColumnFloatMatrix2D((int) v.size(), 1);
V.assign((float[]) v.elements());
} else {
V = factory.appendColumn(V, v);
}
if (C == null) {
C = new DenseFloatMatrix2D(2, 1);
C.assign(alphaBeta);
} else {
C = factory.composeBidiagonal(C, alphaBeta);
}
} else {
FloatMatrix1D u = null;
FloatMatrix1D v = null;
FloatMatrix1D row = null;
if (counter == 1) {
row = U.viewRow(0).copy();
row = M.transApply(row, row);
v = A.zMult(row, v, 1, 0, true);
} else {
row = U.viewRow(0).copy();
row = M.transApply(row, row);
v = A.zMult(row, v, 1, 0, true);
row = V.viewColumn(counter - 2);
v.assign(row, FloatFunctions.plusMultSecond(-C.getQuick(counter - 1, counter - 2)));
}
float alpha = alg.norm2(v);
v.assign(FloatFunctions.div(alpha));
row = A.zMult(v, row);
u = M.apply(row, u);
row = U.viewRow(0);
u.assign(row, FloatFunctions.plusMultSecond(-alpha));
float beta = alg.norm2(u);
alphaBeta.setQuick(0, 0, alpha);
alphaBeta.setQuick(1, 0, beta);
u.assign(FloatFunctions.div(beta));
U.viewRow(0).assign(u);
if (V == null) {
V = new DenseColumnFloatMatrix2D((int) v.size(), 1);
V.assign((float[]) v.elements());
} else {
V = factory.appendColumn(V, v);
}
if (C == null) {
C = new DenseFloatMatrix2D(2, 1);
C.assign(alphaBeta);
} else {
C = factory.composeBidiagonal(C, alphaBeta);
}
counter++;
}
}
public FloatMatrix2D getC() {
return C;
}
public FloatMatrix2D getU() {
return U;
}
public FloatMatrix2D getV() {
return V;
}
}
}
© 2015 - 2025 Weber Informatics LLC | Privacy Policy