package scientific.linear
import std.math.*
import std.unittest.*
import std.unittest.testmacro.*
import scientific.numbers.*
/* Compute the QR factorization of a matrix A. This functions returns
information about elementary reflectors making up of Q.
*/
foreign func LAPACKE_dgeqrf(
matrix_layout: IntNative, // row or column major
m: IntNative, // number of rows in A
n: IntNative, // number of columns in A
a: CPointer<Unit>, // the m-by-n matrix (output - elementary reflectors, Float64)
lda: IntNative, // leading dimension of A
tau: CPointer<Unit> // output - scalar factors of elementary reflectors, Float64
): Int64
/* computes an RQ factorization of a real M-by-N matrix A:
A = R * Q.
*/
foreign func LAPACKE_dgerqf(
matrix_layout: IntNative, // row or column major
m: IntNative, // number of rows in A
n: IntNative, // number of columns in A
a: CPointer<Unit>, // the m-by-n matrix (output - elementary reflectors, Float64)
lda: IntNative, // leading dimension of A
tau: CPointer<Unit> // output - scalar factors of elementary reflectors, Float64
): Int64
/* Applies the result of QR factorization from dgeqrf. */
foreign func LAPACKE_dormqr(
matrix_layout: IntNative, // row or column major
side: UInt8, // side to multiply Q ('l' or 'r')
trans: UInt8, // whether to apply transpose ('n' or 't')
m: IntNative, // number of rows in C
n: IntNative, // number of columns in C
k: IntNative, // number of reflectors in Q
a: CPointer<Unit>, // reflectors as returned from dgeqrf (Float64)
lda: IntNative, // leading dimension of A
tau: CPointer<Unit>, // scalar factors of elementary reflectors (Float64)
c: CPointer<Unit>, // the m-by-n matrix C (output - product of Q with C, Float64)
ldc: IntNative // leading dimension of C
): Int64
/* Applies the result of RQ factorization from dgerqf. */
foreign func LAPACKE_dorgrq(
matrix_layout: IntNative, // row or column major
m: IntNative, // number of rows in C
n: IntNative, // number of columns in C
k: IntNative, // number of reflectors in Q
a: CPointer<Unit>, // reflectors as returned from dgeqrf (Float64)
lda: IntNative, // leading dimension of A
tau: CPointer<Unit> // scalar factors of elementary reflectors (Float64)
): Int64
public func cj_dgeqrf(a: Matrix<Float64>): Vector<Float64> {
let LAPACK_COL_MAJOR: Int64 = 102
let m = a.getRows()
let n = a.getCols()
let lda: Int64 = m
var tau = zeros<Float64>(n)
let info = unsafe {
LAPACKE_dgeqrf(IntNative(LAPACK_COL_MAJOR), IntNative(m), IntNative(n), a.ptr,
IntNative(lda), tau.ptr)
}
return tau
}
public func cj_dgerqf(a: Matrix<Float64>): Vector<Float64> {
let LAPACK_COL_MAJOR: Int64 = 102
let m = a.getRows()
let n = a.getCols()
let minmn = min(m,n)
let lda: Int64 = m
var tau = zeros<Float64>(minmn)
let info = unsafe {
LAPACKE_dgerqf(IntNative(LAPACK_COL_MAJOR), IntNative(m), IntNative(n), a.ptr,
IntNative(lda), tau.ptr)
}
return tau
}
public func rq(a: Matrix<Float64>): (Matrix<Float64>, Matrix<Float64>) {
let r = a.copy()
var tau = cj_dgerqf(r)
let m = a.getRows()
let n = a.getCols()
let minmn = min(m,n)
var q = empty<Float64>(n, n)
if (m >= n) {
for (i in 0..n) {
for (j in 0..n) {
q[i,j] = r[i+m-minmn,j]
}
}
} else {
for (i in 0..m) {
for (j in 0..n) {
q[i+n-minmn,j] = r[i,j]
}
}
}
let info = cj_dorgrq(q, tau)
// change the lower to zero to show R
if (m >= n) {
for (i in 0..m) {
for (j in 0..n) {
if (i > j+m-minmn) {
r[i,j] = 0.0
}
}
}
} else {
for (i in 0..m) {
for (j in 0..n) {
if (i+n-minmn > j) {
r[i,j] = 0.0
}
}
}
}
return (r,q)
}
func testDorgrq(): Unit {
var a = matrix<Float64>(
[[6.0, 2.0, 3.0],
[2.0, 3.0, 1.0]]
)
var tau = cj_dgerqf(a)
var q = empty<Float64>(3, 3)
for (i in 0..2) {
for (j in 0..3) {
q[i+1,j] = a[i,j]
}
}
let LAPACK_COL_MAJOR: Int64 = 102
let m = 3
let n = 3
let lda: Int64 = 3
let k: Int64 = 2
let info = unsafe {
LAPACKE_dorgrq(IntNative(LAPACK_COL_MAJOR), IntNative(m),
IntNative(n), IntNative(k), q.ptr, IntNative(lda), tau.ptr)
}
var expected_q = matrix<Float64>(
[[0.447214, 0.000000, -0.894427],
[0.717137, -0.597614, 0.358569],
[-0.534522, -0.801784, -0.267261]]
)
assertApproxEqual(q, expected_q, atol:1e-6)
}
public func cj_dormqr(a: Matrix<Float64>, tau: Vector<Float64>): Matrix<Float64> {
let LAPACK_COL_MAJOR: Int64 = 102
let side: Rune = r'l'
let trans: Rune = r'n'
let cside: UInt32 = UInt32(side)
let ctrans: UInt32 = UInt32(trans)
let m = a.getRows()
let n = a.getCols()
let lda: Int64 = m
let k: Int64 = m
let ldc: Int64 = m
let c = eye<Float64>(m, n)
let info = unsafe {
LAPACKE_dormqr(IntNative(LAPACK_COL_MAJOR), UInt8(cside), UInt8(ctrans), IntNative(m),
IntNative(n), IntNative(k), a.ptr, IntNative(lda), tau.ptr, c.ptr, IntNative(ldc))
}
return c
}
public func cj_dorgrq(a: Matrix<Float64>, tau: Vector<Float64>) {
let LAPACK_COL_MAJOR: Int64 = 102
let m = a.getRows()
let n = m
let lda: Int64 = m
let k: Int64 = tau.size()
let info = unsafe {
LAPACKE_dorgrq(IntNative(LAPACK_COL_MAJOR), IntNative(m),
IntNative(n), IntNative(k), a.ptr, IntNative(lda), tau.ptr)
}
return info
}
public func cj_dormqr_c(a: Matrix<Float64>, c: Matrix<Float64>, tau: Vector<Float64>): Matrix<Float64> {
let LAPACK_COL_MAJOR: Int64 = 102
let side: Rune = r'l'
let trans: Rune = r'n'
let cside: UInt32 = UInt32(side)
let ctrans: UInt32 = UInt32(trans)
let m = a.getRows()
let n = a.getCols()
let lda: Int64 = m
let k: Int64 = m
let ldc: Int64 = m
let info = unsafe {
LAPACKE_dormqr(IntNative(LAPACK_COL_MAJOR), UInt8(cside), UInt8(ctrans), IntNative(m),
IntNative(n), IntNative(k), a.ptr, IntNative(lda), tau.ptr, c.ptr, IntNative(ldc))
}
return c
}
public func qr(a: Matrix<Float64>): (Matrix<Float64>, Matrix<Float64>) {
let m = a.getRows()
let n = a.getCols()
let acopy = a.copy()
let tau = cj_dgeqrf(acopy)
var res = cj_dormqr(acopy, tau)
// change size of c to m * min(m, n)
var q = res
if (n != min(m, n)) {
q = empty<Float64>(m, min(m, n))
for (i in 0..m) {
for (j in 0..min(m, n)) {
q[i,j] = res[i,j]
}
}
}
// change the lower to zero to show R
let r = empty<Float64>(min(m, n), n)
for (i in 0..min(m, n)) {
for (j in 0..n) {
if (i > j) {
r[i,j] = 0.0
} else {
r[i,j] = acopy[i,j]
}
}
}
return (q, r)
}
public func qrMultiply(a: Matrix<Float64>, c: Matrix<Float64>): (Matrix<Float64>, Matrix<Float64>) {
let m = a.getRows()
let n = a.getCols()
let acopy = a.copy()
let tau = cj_dgeqrf(acopy)
var qc = cj_dormqr_c(acopy, c, tau)
// change the lower to zero to show R
let r = empty<Float64>(min(m, n), n)
for (i in 0..min(m, n)) {
for (j in 0..n) {
if (i > j) {
r[i,j] = 0.0
} else {
r[i,j] = acopy[i,j]
}
}
}
return (qc, r)
}
// Given a and b, return (c, s) for the Givens rotation
public func givens(a: Float64, b: Float64): (Float64, Float64) {
if (b == 0.0) {
return (1.0, 0.0)
} else {
if (abs(b) > abs(a)) {
let tau = -a / b
let s = 1.0 / sqrt(1.0 + tau * tau)
let c = s * tau
return (c, s)
} else {
let tau = -b / a
let c = 1.0 / sqrt(1.0 + tau * tau)
let s = c * tau
return (c, s)
}
}
}
public func givensrotation(q1: Matrix<Float64>,r1: Matrix<Float64>,h: Int64, l: Int64): (Matrix<Float64>, Matrix<Float64>) {
let n = r1.getCols()
let m1 = q1.getRows()
let (c, s) = givens(r1[h-1,l], r1[h,l])
for (j in 0..m1){
let tau1 = q1[j,l]
let tau2 = q1[j,h]
q1[j,l] = c*tau1 - s*tau2
q1[j,h] = s*tau1 + c*tau2
}
for (j in 0..n){
let tau1 = r1[l,j]
let tau2 = r1[h,j]
r1[l,j] = c*tau1 - s*tau2
r1[h,j] = s*tau1 + c*tau2
}
return (q1, r1)
}
public func givensrotation_qn(qn: Matrix<Float64>, q1: Matrix<Float64>, r1: Matrix<Float64>, h: Int64, l: Int64): (Matrix<Float64>, Matrix<Float64>, Matrix<Float64>) {
let (c, s) = givens(qn[h-1,l], qn[h,l])
for(j in 0..qn.getCols()) {
let tau1 = qn[h-1,j]
let tau2 = qn[h,j]
qn[h-1,j] = c*tau1 - s*tau2
qn[h,j] = s*tau1 + c*tau2
}
for(j in 0..r1.getCols()) {
let tau1 = r1[h-1,j]
let tau2 = r1[h,j]
r1[h-1,j] = c*tau1 - s*tau2
r1[h,j] = s*tau1 + c*tau2
}
for(j in 0..q1.getRows()) {
let tau1 = q1[j,h-1]
let tau2 = q1[j,h]
q1[j,h-1] = c*tau1 - s*tau2
q1[j,h] = s*tau1 + c*tau2
}
return (qn, q1, r1)
}
public func qrDelete(q: Matrix<Float64>, r: Matrix<Float64>, k: Int64, p: Int64, alter: String): (Matrix<Float64>, Matrix<Float64>) {
let m1 = q.getRows()
let n1 = q.getCols()
var r2 = empty<Float64>(r.getRows(), r.getCols())
var qn = empty<Float64>(q.getCols(), p)
var q2 = q.copy()
let m2 = r.getRows()
let n2 = r.getCols()
if(alter == "col"){
r2 = empty<Float64>(r.getRows(), r.getCols() - p)
for (i in 0..r.getRows()) {
for (j in 0..k) {
r2[i,j] = r[i,j]
}
for (j in k+p..r.getCols()) {
r2[i,j-p] = r[i,j]
}
}
q2 = q.copy()
let m = r2.getRows()
let n = r2.getCols()
for (j in 0..n) {
for (i in m-1..0:-1) {
if(i>j && r2[i,j]!=0.000000){
var (q2new, r2new) = givensrotation(q2, r2, i, j)
q2 = q2new
r2 = r2new
}
}
}
}else if(alter == "row"){
var qn = empty<Float64>(q.getCols(), p)
var q1 = q.copy()
for(i in 0..p) {
for(j in 0..q.getCols()) {
qn[j,i] = q[i+k,j]
}
}
var r1 = r.copy()
for(i in 0..qn.getCols()) {
for(j in qn.getRows()-1..0:-1) {
var (qnnew, q1new, r1new) = givensrotation_qn(qn, q1, r1, j, i)
qn = qnnew
q1 = q1new
r1 = r1new
}
}
q2 = empty<Float64>(q.getRows()-p, q.getCols()-p)
for(i in 0..k) {
for(j in p..q1.getCols()) {
q2[i,j-p] = q1[i,j]
}
}
for(i in k+p..q1.getRows()) {
for(j in p..q1.getCols()) {
q2[i-p,j-p] = q1[i,j]
}
}
r2 = empty<Float64>(r.getRows()-p, r.getCols())
for(i in 0..r1.getCols()) {
for(j in p..r1.getRows()) {
r2[j-p,i] = r1[j,i]
}
}
}
return (q2,r2)
}
public func qrInsert(q: Matrix<Float64>, r: Matrix<Float64>, u: Matrix<Float64>, k: Int64, alter: String): (Matrix<Float64>, Matrix<Float64>) {
let m1 = q.getRows()
let n1 = q.getCols()
var r2 = empty<Float64>(r.getRows(), r.getCols())
var q2 = q.copy()
let m2 = r.getRows()
let n2 = r.getCols()
println(q)
println(r)
var qt = q.transpose()
println(qt)
var w = qt*u
println(w)
if(alter == "col"){
r2 = empty<Float64>(r.getRows(), r.getCols() + w.getCols())
for (i in 0..r.getRows()) {
for (j in 0..k) {
r2[i,j] = r[i,j]
}
for (j in k..k+w.getCols()) {
r2[i,j] = w[i,j-k]
}
for(j in k..r.getCols()) {
r2[i,j+w.getCols()] = r[i,j]
}
}
println(r2)
q2 = q.copy()
let m = r2.getRows()
let n = r2.getCols()
for (j in 0..n) {
for (i in m-1..0:-1) {
if(i>j && r2[i,j]!=0.000000){
var (q2new, r2new) = givensrotation(q2, r2, i, j)
q2 = q2new
r2 = r2new
}
}
}
}else if(alter == "row"){
}
println(q2)
println(r2)
return (q2,r2)
}
func testDormqr(): Unit {
var a = matrix<Float64>(
[[1.0, 0.0, 0.0],
[1.0, 1.0, 1.0],
[0.0, 0.0, 1.0]]
)
let m = a.getRows()
let n = a.getCols()
let lda: Int64 = m
var tau = zeros<Float64>(n)
var cinfo = unsafe { LAPACKE_dgeqrf(LAPACK_COL_MAJOR, IntNative(m), IntNative(n),
a.ptr, IntNative(lda), tau.ptr) }
let side: Rune = r'l'
let trans: Rune = r'n'
let cside: UInt32 = UInt32(side)
let ctrans: UInt32 = UInt32(trans)
let k: Int64 = m
let ldc: Int64 = m
var c = matrix<Float64>(
[[1.0, 0.0, 0.0],
[0.0, 1.0, 0.0],
[0.0, 0.0, 1.0]]
)
var info = unsafe { LAPACKE_dormqr(LAPACK_COL_MAJOR, UInt8(cside), UInt8(ctrans),
IntNative(m), IntNative(n), IntNative(k), a.ptr, IntNative(lda),
tau.ptr, c.ptr, IntNative(ldc)) }
// change the lower to zero to show R
for (i in 0..m) {
for (j in 0..n) {
if (i > j) {
a[i,j] = 0.0
}
}
}
var expected_factor = matrix<Float64>(
[[-1.414214, -0.707107, -0.707107],
[0.000000, 0.707107, 0.707107],
[0.000000, 0.000000, 1.000000]]
)
@Assert(approxEqual(a, expected_factor, atol:1e-6))
}
@Test
public class TestQR {
@TestCase
func testDgeqrf(): Unit {
var a = matrix<Float64>(
[[0.0, 3.0, 1.0],
[0.0, 4.0, -2.0],
[2.0, 1.0, 2.0]]
)
let m = a.getRows()
let n = a.getCols()
var info = cj_dgeqrf(a)
var expected_factor = matrix<Float64>(
[[-2.000000, -1.000000, -2.000000],
[0.000000, -5.000000, 1.000000],
[1.000000, -0.333333, -2.000000]]
)
@Assert(approxEqual(a, expected_factor, atol:1e-6))
}
@TestCase
func testLapackDormqr(): Unit {
testDormqr()
}
@TestCase
func testQr1(): Unit {
var a = matrix<Float64>(
[[1.0, 0.0, 0.0],
[1.0, 1.0, 1.0],
[0.0, 0.0, 1.0]]
)
let (q, r) = qr(a)
var expected_q = matrix<Float64>(
[[-0.707107, -0.707107, 0.000000],
[-0.707107, 0.707107, 0.000000],
[ 0.000000, 0.000000, 1.000000]]
)
var expected_r = matrix<Float64>(
[[-1.414214, -0.707107, -0.707107],
[ 0.000000, 0.707107, 0.707107],
[ 0.000000, 0.000000, 1.000000]]
)
@Assert(approxEqual(q, expected_q, atol:1e-6))
@Assert(approxEqual(r, expected_r, atol:1e-6))
@Assert(approxEqual(q * r, a, atol:1e-6))
}
@TestCase
func testQr2(): Unit {
var a = matrix<Float64>(
[[1.0, 3.0, 3.0],
[2.0, 3.0, 2.0],
[2.0, 3.0, 3.0],
[1.0, 3.0, 2.0]]
)
let (q, r) = qr(a)
var expected_q = matrix<Float64>(
[[-0.316228, 0.632456, -0.500000],
[-0.632456, -0.316228, 0.500000],
[-0.632456, -0.316228, -0.500000],
[-0.316228, 0.632456, 0.500000]]
)
var expected_r = matrix<Float64>(
[[-3.162278, -5.692100, -4.743416],
[0.000000, 1.897367, 1.581139],
[0.000000, 0.000000, -1.000000]]
)
@Assert(approxEqual(q, expected_q, atol:1e-6))
@Assert(approxEqual(r, expected_r, atol:1e-6))
@Assert(approxEqual(q * r, a, atol:1e-6))
}
@TestCase
func testQr3(): Unit {
var a = matrix<Float64>(
[[1.0, 3.0, 3.0],
[2.0, 3.0, 2.0],
[2.0, 3.0, 3.0],
[1.0, 3.0, 2.0]]
).transpose()
let (q, r) = qr(a)
var expected_q = matrix<Float64>(
[[-0.229416, 0.826234, 0.514496],
[-0.688247, 0.236067, -0.685994],
[-0.688247, -0.511478, 0.514496]]
)
var expected_r = matrix<Float64>(
[[-4.358899, -3.900067, -4.588315, -3.670652],
[0.000000, 1.337712, 0.826234, 0.511478],
[0.000000, 0.000000, 0.514496, -0.514496]]
)
@Assert(approxEqual(q, expected_q, atol:1e-6))
@Assert(approxEqual(r, expected_r, atol:1e-6))
@Assert(approxEqual(q * r, a, atol:1e-6))
}
@TestCase
func testQrMultiply(): Unit {
var a = matrix<Float64>(
[[1.0, 3.0, 3.0],
[2.0, 3.0, 2.0],
[2.0, 3.0, 3.0],
[1.0, 3.0, 2.0]]
)
let (qc, r) = qrMultiply(a, 2.0 * eye(4, 3))
var expected_qc = matrix<Float64>(
[[-0.632456, 1.264911, -1.000000],
[-1.264911, -0.632456, 1.000000],
[-1.264911, -0.632456, -1.000000],
[-0.632456, 1.264911, 1.000000]]
)
var expected_r = matrix<Float64>(
[[-3.162278, -5.692100, -4.743416],
[0.000000, 1.897367, 1.581139],
[0.000000, 0.000000, -1.000000]]
)
@Assert(approxEqual(qc, expected_qc, atol:1e-6))
@Assert(approxEqual(r, expected_r, atol:1e-6))
}
@TestCase
func testQrDelete1(): Unit {
var a = matrix<Float64>(
[[ 3.0, -2.0, -2.0, 1.0],
[ 6.0, -9.0, -3.0, 2.0],
[ -3.0, 10.0, 1.0, 3.0]])
let (q, r) = qr(a)
let (q1, r1) = qrDelete(q, r, 1, 2, "col")
var expected_q = matrix<Float64>(
[[-0.408248, 0.182574, -0.894427],
[-0.816497, 0.365148, 0.447214],
[0.408248, 0.912871, 0.000000]]
)
var expected_r = matrix<Float64>(
[[-7.348469, -0.816497],
[0.000000, 3.651484],
[0.000000, 0.000000]]
)
@Assert(approxEqual(q1, expected_q, atol:1e-6))
@Assert(approxEqual(r1, expected_r, atol:1e-6))
}
@TestCase
func testQrDelete2(): Unit {
var a = matrix<Float64>(
[[ 3.0, -2.0, -2.0, 1.0],
[ 6.0, -9.0, -3.0, 2.0],
[ -3.0, 10.0, 1.0, 3.0]])
let (q, r) = qr(a)
let (q1, r1) = qrDelete(q, r, 0, 1, "col")
var expected_q = matrix<Float64>(
[[-0.147043, 0.702303, -0.696526],
[-0.661693, 0.453571, 0.597022],
[0.735215, 0.548674, 0.398015]]
)
var expected_r = matrix<Float64>(
[[13.601471, 3.014380, 0.735215],
[0.000000, -2.216645, 3.255468],
[0.000000, 0.000000, 1.691563]]
)
@Assert(approxEqual(q1, expected_q, atol:1e-6))
@Assert(approxEqual(r1, expected_r, atol:1e-6))
}
@TestCase
func testQrDelete3(): Unit {
var a = matrix<Float64>(
[[ 3.0, -2.0, -2.0, 1.0],
[ 6.0, -9.0, -3.0, 2.0],
[ -3.0, 10.0, 1.0, 3.0]]
)
let (q, r) = qr(a)
let (q1, r1) = qrDelete(q, r, 0, 1, "row")
var expected_q = matrix<Float64>(
[[-0.894427, 0.447214],
[0.447214, 0.894427]]
)
var expected_r = matrix<Float64>(
[[-6.708204, 12.521981, 3.130495, -0.447214],
[0.000000, 4.919350, -0.447214, 3.577709]]
)
@Assert(approxEqual(q1, expected_q, atol:1e-6))
@Assert(approxEqual(r1, expected_r, atol:1e-6))
}
@TestCase
func testQrDelete4(): Unit {
var a = matrix<Float64>(
[[ 3.0, -2.0, -2.0, 1.0],
[ 6.0, -9.0, -3.0, 2.0],
[ -3.0, 10.0, 1.0, 3.0]]
)
let (q, r) = qr(a)
let (q1, r1) = qrDelete(q, r, 1, 2, "row")
var expected_q = matrix<Float64>(
[[1.000000]]
)
var expected_r = matrix<Float64>(
[[3.000000, -2.000000, -2.000000, 1.000000]]
)
@Assert(approxEqual(q1, expected_q, atol:1e-6))
@Assert(approxEqual(r1, expected_r, atol:1e-6))
}
@TestCase
func testDgerqf(): Unit {
var a = matrix<Float64>(
[[6.0, 2.0, 3.0],
[2.0, 3.0, 1.0]]
)
var tau = cj_dgerqf(a)
var expected_tau = vector<Float64>([1.824477, 1.267261])
@Assert(approxEqual(tau, expected_tau, atol:1e-6))
}
@TestCase
func testLapackDorgrq(): Unit {
testDorgrq()
}
@TestCase
func testrq(): Unit {
var a = matrix<Float64>(
[[6.0, 2.0, 3.0],
[2.0, 3.0, 1.0]]
)
let (r, q) = rq(a)
var expected_r = matrix<Float64>(
[[0.000000, 4.183300, -5.612486],
[0.000000, 0.000000, -3.741657]]
)
var expected_q = matrix<Float64>(
[[0.447214, 0.000000, -0.894427],
[0.717137, -0.597614, 0.358569],
[-0.534522, -0.801784, -0.267261]]
)
@Assert(approxEqual(r, expected_r, atol:1e-6))
@Assert(approxEqual(q, expected_q, atol:1e-6))
}
@TestCase
func testrq2(): Unit {
var a = matrix<Float64>(
[[6.0, 2.0, 3.0],
[2.0, 3.0, 1.0],
[3.0, 4.0, 8.0],
[4.0, 7.0, 9.0]]
)
let (r, q) = rq(a)
var expected_r = matrix<Float64>(
[[4.478343, 0.078027, -5.379438],
[0.942809, -1.794631, -3.144902],
[0.000000, 1.755617, -9.269186],
[0.000000, 0.000000, -12.083046]]
)
var expected_q = matrix<Float64>(
[[0.942809, -0.235702, -0.235702],
[-0.039013, -0.780274, 0.624219],
[-0.331042, -0.579324, -0.744845]]
)
@Assert(approxEqual(r, expected_r, atol:1e-6))
@Assert(approxEqual(q, expected_q, atol:1e-6))
}
@TestCase
func testQrInsert1(): Unit {
var a = matrix<Float64>(
[[ 3.0, -2.0, -2.0],
[ 6.0, -9.0, -3.0],
[ -3.0, 10.0, 1.0]]
)
var u = matrix<Float64>(
[[ 6.0],
[ -3.0],
[ 9.0]]
)
let (q, r) = qr(a)
let (q1, r1) = qrInsert(q, r, u, 1, "col")
var expected_q = matrix<Float64>(
[[-0.408248, 0.707107, -0.577350],
[-0.816497, 0.000000, 0.577350],
[0.408248, 0.707107, 0.577350]]
)
var expected_r = matrix<Float64>(
[[-7.348469, 3.674235, 12.247449, 3.674235],
[0.000000, 10.606602, 5.656854, -0.707107],
[0.000000, 0.000000, 1.732051, -0.000000]]
)
@Assert(approxEqual(q1, expected_q, atol:1e-6))
@Assert(approxEqual(r1, expected_r, atol:1e-6))
}
@TestCase
func testQrInsert2(): Unit {
var a = matrix<Float64>(
[[ 3.0, -2.0, -2.0],
[ 6.0, -9.0, -3.0],
[ -3.0, 10.0, 1.0]]
)
var u = matrix<Float64>(
[[ 1.0, 3.0],
[ 2.0, -2.0],
[ 3.0, 1.0]]
)
let (q, r) = qr(a)
let (q1, r1) = qrInsert(q, r, u, 2, "col")
var expected_q = matrix<Float64>(
[[-0.408248, 0.507093, -0.759072],
[-0.816497, 0.169031, 0.552052],
[0.408248, 0.845154, 0.345033]]
)
var expected_r = matrix<Float64>(
[[-7.348469, 12.247449, -0.816497, 0.816497, 3.674235],
[0.000000, 5.916080, 3.380617, 2.028370, -0.676123],
[0.000000, 0.000000, 1.380131, -3.036288, 0.207020]]
)
@Assert(approxEqual(q1, expected_q, atol:1e-6))
@Assert(approxEqual(r1, expected_r, atol:1e-6))
}
}