mirror of
				https://github.com/gonum/gonum.git
				synced 2025-10-31 18:42:45 +08:00 
			
		
		
		
	matrix/mat64: reverse API signatures for matrix extraction
This commit is contained in:
		| @@ -203,9 +203,9 @@ func (gsvd *GSVD) ValuesB(s []float64) []float64 { | ||||
| 	return s | ||||
| } | ||||
|  | ||||
| // ZeroRFromGSVD extracts the matrix [ 0 R ] from the singular value decomposition, storing | ||||
| // the result in-place into the receiver. [ 0 R ] is size (k+l)×c. | ||||
| func (m *Dense) ZeroRFromGSVD(gsvd *GSVD) { | ||||
| // ZeroRTo extracts the matrix [ 0 R ] from the singular value decomposition, storing | ||||
| // the result in-place into dst. [ 0 R ] is size (k+l)×c. | ||||
| func (gsvd *GSVD) ZeroRTo(dst *Dense) { | ||||
| 	if gsvd.kind == 0 { | ||||
| 		panic("gsvd: no decomposition computed") | ||||
| 	} | ||||
| @@ -214,13 +214,13 @@ func (m *Dense) ZeroRFromGSVD(gsvd *GSVD) { | ||||
| 	k := gsvd.k | ||||
| 	l := gsvd.l | ||||
| 	h := min(k+l, r) | ||||
| 	m.reuseAsZeroed(k+l, c) | ||||
| 	dst.reuseAsZeroed(k+l, c) | ||||
| 	a := Dense{ | ||||
| 		mat:     gsvd.a, | ||||
| 		capRows: r, | ||||
| 		capCols: c, | ||||
| 	} | ||||
| 	m.Slice(0, h, c-k-l, c).(*Dense). | ||||
| 	dst.Slice(0, h, c-k-l, c).(*Dense). | ||||
| 		Copy(a.Slice(0, h, c-k-l, c)) | ||||
| 	if r < k+l { | ||||
| 		b := Dense{ | ||||
| @@ -228,32 +228,32 @@ func (m *Dense) ZeroRFromGSVD(gsvd *GSVD) { | ||||
| 			capRows: gsvd.p, | ||||
| 			capCols: c, | ||||
| 		} | ||||
| 		m.Slice(r, k+l, c+r-k-l, c).(*Dense). | ||||
| 		dst.Slice(r, k+l, c+r-k-l, c).(*Dense). | ||||
| 			Copy(b.Slice(r-k, l, c+r-k-l, c)) | ||||
| 	} | ||||
| } | ||||
|  | ||||
| // SigmaAFromGSVD extracts the matrix Σ₁ from the singular value decomposition, storing | ||||
| // the result in-place into the receiver. Σ₁ is size r×(k+l). | ||||
| func (m *Dense) SigmaAFromGSVD(gsvd *GSVD) { | ||||
| // SigmaATo extracts the matrix Σ₁ from the singular value decomposition, storing | ||||
| // the result in-place into dst. Σ₁ is size r×(k+l). | ||||
| func (gsvd *GSVD) SigmaATo(dst *Dense) { | ||||
| 	if gsvd.kind == 0 { | ||||
| 		panic("gsvd: no decomposition computed") | ||||
| 	} | ||||
| 	r := gsvd.r | ||||
| 	k := gsvd.k | ||||
| 	l := gsvd.l | ||||
| 	m.reuseAsZeroed(r, k+l) | ||||
| 	dst.reuseAsZeroed(r, k+l) | ||||
| 	for i := 0; i < k; i++ { | ||||
| 		m.set(i, i, 1) | ||||
| 		dst.set(i, i, 1) | ||||
| 	} | ||||
| 	for i := k; i < min(r, k+l); i++ { | ||||
| 		m.set(i, i, gsvd.s1[i]) | ||||
| 		dst.set(i, i, gsvd.s1[i]) | ||||
| 	} | ||||
| } | ||||
|  | ||||
| // SigmaBFromGSVD extracts the matrix Σ₂ from the singular value decomposition, storing | ||||
| // the result in-place into the receiver. Σ₂ is size p×(k+l). | ||||
| func (m *Dense) SigmaBFromGSVD(gsvd *GSVD) { | ||||
| // SigmaBTo extracts the matrix Σ₂ from the singular value decomposition, storing | ||||
| // the result in-place into dst. Σ₂ is size p×(k+l). | ||||
| func (gsvd *GSVD) SigmaBTo(dst *Dense) { | ||||
| 	if gsvd.kind == 0 { | ||||
| 		panic("gsvd: no decomposition computed") | ||||
| 	} | ||||
| @@ -261,65 +261,65 @@ func (m *Dense) SigmaBFromGSVD(gsvd *GSVD) { | ||||
| 	p := gsvd.p | ||||
| 	k := gsvd.k | ||||
| 	l := gsvd.l | ||||
| 	m.reuseAsZeroed(p, k+l) | ||||
| 	dst.reuseAsZeroed(p, k+l) | ||||
| 	for i := 0; i < min(l, r-k); i++ { | ||||
| 		m.set(i, i+k, gsvd.s2[k+i]) | ||||
| 		dst.set(i, i+k, gsvd.s2[k+i]) | ||||
| 	} | ||||
| 	for i := r - k; i < l; i++ { | ||||
| 		m.set(i, i+k, 1) | ||||
| 		dst.set(i, i+k, 1) | ||||
| 	} | ||||
| } | ||||
|  | ||||
| // UFromGSVD extracts the matrix U from the singular value decomposition, storing | ||||
| // the result in-place into the receiver. U is size r×r. | ||||
| func (m *Dense) UFromGSVD(gsvd *GSVD) { | ||||
| // UTo extracts the matrix U from the singular value decomposition, storing | ||||
| // the result in-place into dst. U is size r×r. | ||||
| func (gsvd *GSVD) UTo(dst *Dense) { | ||||
| 	if gsvd.kind&matrix.GSVDU == 0 { | ||||
| 		panic("mat64: improper GSVD kind") | ||||
| 	} | ||||
| 	r := gsvd.u.Rows | ||||
| 	c := gsvd.u.Cols | ||||
| 	m.reuseAs(r, c) | ||||
| 	dst.reuseAs(r, c) | ||||
|  | ||||
| 	tmp := &Dense{ | ||||
| 		mat:     gsvd.u, | ||||
| 		capRows: r, | ||||
| 		capCols: c, | ||||
| 	} | ||||
| 	m.Copy(tmp) | ||||
| 	dst.Copy(tmp) | ||||
| } | ||||
|  | ||||
| // VFromGSVD extracts the matrix V from the singular value decomposition, storing | ||||
| // the result in-place into the receiver. V is size p×p. | ||||
| func (m *Dense) VFromGSVD(gsvd *GSVD) { | ||||
| // VTo extracts the matrix V from the singular value decomposition, storing | ||||
| // the result in-place into dst. V is size p×p. | ||||
| func (gsvd *GSVD) VTo(dst *Dense) { | ||||
| 	if gsvd.kind&matrix.GSVDV == 0 { | ||||
| 		panic("mat64: improper GSVD kind") | ||||
| 	} | ||||
| 	r := gsvd.v.Rows | ||||
| 	c := gsvd.v.Cols | ||||
| 	m.reuseAs(r, c) | ||||
| 	dst.reuseAs(r, c) | ||||
|  | ||||
| 	tmp := &Dense{ | ||||
| 		mat:     gsvd.v, | ||||
| 		capRows: r, | ||||
| 		capCols: c, | ||||
| 	} | ||||
| 	m.Copy(tmp) | ||||
| 	dst.Copy(tmp) | ||||
| } | ||||
|  | ||||
| // QFromGSVD extracts the matrix Q from the singular value decomposition, storing | ||||
| // the result in-place into the receiver. Q is size c×c. | ||||
| func (m *Dense) QFromGSVD(gsvd *GSVD) { | ||||
| // QTo extracts the matrix Q from the singular value decomposition, storing | ||||
| // the result in-place into dst. Q is size c×c. | ||||
| func (gsvd *GSVD) QTo(dst *Dense) { | ||||
| 	if gsvd.kind&matrix.GSVDQ == 0 { | ||||
| 		panic("mat64: improper GSVD kind") | ||||
| 	} | ||||
| 	r := gsvd.q.Rows | ||||
| 	c := gsvd.q.Cols | ||||
| 	m.reuseAs(r, c) | ||||
| 	dst.reuseAs(r, c) | ||||
|  | ||||
| 	tmp := &Dense{ | ||||
| 		mat:     gsvd.q, | ||||
| 		capRows: r, | ||||
| 		capCols: c, | ||||
| 	} | ||||
| 	m.Copy(tmp) | ||||
| 	dst.Copy(tmp) | ||||
| } | ||||
|   | ||||
| @@ -26,8 +26,8 @@ func ExampleGSVD() { | ||||
| 	} | ||||
|  | ||||
| 	var u, v mat64.Dense | ||||
| 	u.UFromGSVD(&gsvd) | ||||
| 	v.VFromGSVD(&gsvd) | ||||
| 	gsvd.UTo(&u) | ||||
| 	gsvd.VTo(&v) | ||||
|  | ||||
| 	s1 := gsvd.ValuesA(nil) | ||||
| 	s2 := gsvd.ValuesB(nil) | ||||
| @@ -38,8 +38,8 @@ func ExampleGSVD() { | ||||
| 		s2, mat64.Formatted(&v, mat64.Prefix("\t    "), mat64.Excerpt(2))) | ||||
|  | ||||
| 	var zeroR, q mat64.Dense | ||||
| 	zeroR.ZeroRFromGSVD(&gsvd) | ||||
| 	q.QFromGSVD(&gsvd) | ||||
| 	gsvd.ZeroRTo(&zeroR) | ||||
| 	gsvd.QTo(&q) | ||||
| 	q.Mul(&zeroR, &q) | ||||
| 	fmt.Printf("\nCommon basis vectors\n\n\tQ^T = %.4f\n", | ||||
| 		mat64.Formatted(q.T(), mat64.Prefix("\t      "))) | ||||
|   | ||||
| @@ -107,12 +107,12 @@ func TestGSVD(t *testing.T) { | ||||
|  | ||||
| func extractGSVD(gsvd *GSVD) (c, s []float64, s1, s2, zR, u, v, q *Dense) { | ||||
| 	var s1m, s2m, zeroR, um, vm, qm Dense | ||||
| 	s1m.SigmaAFromGSVD(gsvd) | ||||
| 	s2m.SigmaBFromGSVD(gsvd) | ||||
| 	zeroR.ZeroRFromGSVD(gsvd) | ||||
| 	um.UFromGSVD(gsvd) | ||||
| 	vm.VFromGSVD(gsvd) | ||||
| 	qm.QFromGSVD(gsvd) | ||||
| 	gsvd.SigmaATo(&s1m) | ||||
| 	gsvd.SigmaBTo(&s2m) | ||||
| 	gsvd.ZeroRTo(&zeroR) | ||||
| 	gsvd.UTo(&um) | ||||
| 	gsvd.VTo(&vm) | ||||
| 	gsvd.QTo(&qm) | ||||
| 	c = gsvd.ValuesA(nil) | ||||
| 	s = gsvd.ValuesB(nil) | ||||
| 	return c, s, &s1m, &s2m, &zeroR, &um, &vm, &qm | ||||
|   | ||||
| @@ -141,11 +141,11 @@ func (gsvd *HOGSVD) Len() int { | ||||
| 	return gsvd.n | ||||
| } | ||||
|  | ||||
| // UFromHOGSVD extracts the matrix U_n from the singular value decomposition, storing | ||||
| // the result in-place into the receiver. U_n is size r×c. | ||||
| // UTo extracts the matrix U_n from the singular value decomposition, storing | ||||
| // the result in-place into dst. U_n is size r×c. | ||||
| // | ||||
| // UFromHOGSVD will panic if the receiver does not contain a successful factorization. | ||||
| func (m *Dense) UFromHOGSVD(gsvd *HOGSVD, n int) { | ||||
| // UTo will panic if the receiver does not contain a successful factorization. | ||||
| func (gsvd *HOGSVD) UTo(dst *Dense, n int) { | ||||
| 	if gsvd.n == 0 { | ||||
| 		panic("hogsvd: unsuccessful factorization") | ||||
| 	} | ||||
| @@ -153,10 +153,10 @@ func (m *Dense) UFromHOGSVD(gsvd *HOGSVD, n int) { | ||||
| 		panic("hogsvd: invalid index") | ||||
| 	} | ||||
|  | ||||
| 	m.reuseAs(gsvd.b[n].Dims()) | ||||
| 	m.Copy(&gsvd.b[n]) | ||||
| 	dst.reuseAs(gsvd.b[n].Dims()) | ||||
| 	dst.Copy(&gsvd.b[n]) | ||||
| 	for j, f := range gsvd.Values(nil, n) { | ||||
| 		v := m.ColView(j) | ||||
| 		v := dst.ColView(j) | ||||
| 		v.ScaleVec(1/f, v) | ||||
| 	} | ||||
| } | ||||
| @@ -188,13 +188,14 @@ func (gsvd *HOGSVD) Values(s []float64, n int) []float64 { | ||||
| 	return s | ||||
| } | ||||
|  | ||||
| // VFromHOGSVD extracts the matrix V from the singular value decomposition, storing | ||||
| // the result in-place into the receiver. V is size c×c. | ||||
| // VTo extracts the matrix V from the singular value decomposition, storing | ||||
| // the result in-place into dst. V is size c×c. | ||||
| // | ||||
| // VFromHOGSVD will panic if the receiver does not contain a successful factorization. | ||||
| func (m *Dense) VFromHOGSVD(gsvd *HOGSVD) { | ||||
| // VTo will panic if the receiver does not contain a successful factorization. | ||||
| func (gsvd *HOGSVD) VTo(dst *Dense) { | ||||
| 	if gsvd.n == 0 { | ||||
| 		panic("hogsvd: unsuccessful factorization") | ||||
| 	} | ||||
| 	*m = *DenseCopyOf(gsvd.v) | ||||
| 	dst.reuseAs(gsvd.v.Dims()) | ||||
| 	dst.Copy(gsvd.v) | ||||
| } | ||||
|   | ||||
| @@ -25,14 +25,14 @@ func ExampleHOGSVD() { | ||||
|  | ||||
| 	for i, n := range []string{"Africa", "Asia", "Latin America/Caribbean", "Oceania"} { | ||||
| 		var u mat64.Dense | ||||
| 		u.UFromHOGSVD(&gsvd, i) | ||||
| 		gsvd.UTo(&u, i) | ||||
| 		s := gsvd.Values(nil, i) | ||||
| 		fmt.Printf("%s\n\ts_%d = %.4f\n\n\tU_%[2]d = %.4[4]f\n", | ||||
| 			n, i, s, mat64.Formatted(&u, mat64.Prefix("\t      "))) | ||||
| 	} | ||||
|  | ||||
| 	var v mat64.Dense | ||||
| 	v.VFromHOGSVD(&gsvd) | ||||
| 	gsvd.VTo(&v) | ||||
| 	fmt.Printf("\nCommon basis vectors\n\n\tV^T = %.4f", | ||||
| 		mat64.Formatted(v.T(), mat64.Prefix("\t      "))) | ||||
|  | ||||
|   | ||||
| @@ -80,10 +80,10 @@ func extractHOGSVD(gsvd *HOGSVD) (u []*Dense, s [][]float64, v *Dense) { | ||||
| 	s = make([][]float64, gsvd.Len()) | ||||
| 	for i := 0; i < gsvd.Len(); i++ { | ||||
| 		u[i] = &Dense{} | ||||
| 		u[i].UFromHOGSVD(gsvd, i) | ||||
| 		gsvd.UTo(u[i], i) | ||||
| 		s[i] = gsvd.Values(nil, i) | ||||
| 	} | ||||
| 	v = &Dense{} | ||||
| 	v.VFromHOGSVD(gsvd) | ||||
| 	gsvd.VTo(v) | ||||
| 	return u, s, v | ||||
| } | ||||
|   | ||||
| @@ -36,7 +36,7 @@ func (lq *LQ) updateCond() { | ||||
| // | ||||
| // The LQ decomposition is a factorization of the matrix A such that A = L * Q. | ||||
| // The matrix Q is an orthonormal n×n matrix, and L is an m×n upper triangular matrix. | ||||
| // L and Q can be extracted from the LFromLQ and QFromLQ methods on Dense. | ||||
| // L and Q can be extracted from the LTo and QTo methods. | ||||
| func (lq *LQ) Factorize(a Matrix) { | ||||
| 	m, n := a.Dims() | ||||
| 	if m > n { | ||||
| @@ -58,10 +58,10 @@ func (lq *LQ) Factorize(a Matrix) { | ||||
| // TODO(btracey): Add in the "Reduced" forms for extracting the m×m orthogonal | ||||
| // and upper triangular matrices. | ||||
|  | ||||
| // LFromLQ extracts the m×n lower trapezoidal matrix from a LQ decomposition. | ||||
| func (m *Dense) LFromLQ(lq *LQ) { | ||||
| // LTo extracts the m×n lower trapezoidal matrix from a LQ decomposition. | ||||
| func (lq *LQ) LTo(dst *Dense) { | ||||
| 	r, c := lq.lq.Dims() | ||||
| 	m.reuseAs(r, c) | ||||
| 	dst.reuseAs(r, c) | ||||
|  | ||||
| 	// Disguise the LQ as a lower triangular | ||||
| 	t := &TriDense{ | ||||
| @@ -74,25 +74,25 @@ func (m *Dense) LFromLQ(lq *LQ) { | ||||
| 		}, | ||||
| 		cap: lq.lq.capCols, | ||||
| 	} | ||||
| 	m.Copy(t) | ||||
| 	dst.Copy(t) | ||||
|  | ||||
| 	if r == c { | ||||
| 		return | ||||
| 	} | ||||
| 	// Zero right of the triangular. | ||||
| 	for i := 0; i < r; i++ { | ||||
| 		zero(m.mat.Data[i*m.mat.Stride+r : i*m.mat.Stride+c]) | ||||
| 		zero(dst.mat.Data[i*dst.mat.Stride+r : i*dst.mat.Stride+c]) | ||||
| 	} | ||||
| } | ||||
|  | ||||
| // QFromLQ extracts the n×n orthonormal matrix Q from an LQ decomposition. | ||||
| func (m *Dense) QFromLQ(lq *LQ) { | ||||
| // QTo extracts the n×n orthonormal matrix Q from an LQ decomposition. | ||||
| func (lq *LQ) QTo(dst *Dense) { | ||||
| 	r, c := lq.lq.Dims() | ||||
| 	m.reuseAs(c, c) | ||||
| 	dst.reuseAs(c, c) | ||||
|  | ||||
| 	// Set Q = I. | ||||
| 	for i := 0; i < c; i++ { | ||||
| 		v := m.mat.Data[i*m.mat.Stride : i*m.mat.Stride+c] | ||||
| 		v := dst.mat.Data[i*dst.mat.Stride : i*dst.mat.Stride+c] | ||||
| 		zero(v) | ||||
| 		v[i] = 1 | ||||
| 	} | ||||
| @@ -127,10 +127,10 @@ func (m *Dense) QFromLQ(lq *LQ) { | ||||
|  | ||||
| 		// Compute the multiplication matrix. | ||||
| 		blas64.Ger(-lq.tau[i], v, v, h) | ||||
| 		qCopy.Copy(m) | ||||
| 		qCopy.Copy(dst) | ||||
| 		blas64.Gemm(blas.NoTrans, blas.NoTrans, | ||||
| 			1, h, qCopy.mat, | ||||
| 			0, m.mat) | ||||
| 			0, dst.mat) | ||||
| 	} | ||||
| } | ||||
|  | ||||
|   | ||||
| @@ -27,16 +27,16 @@ func TestLQ(t *testing.T) { | ||||
| 		var want Dense | ||||
| 		want.Clone(a) | ||||
|  | ||||
| 		lq := &LQ{} | ||||
| 		var lq LQ | ||||
| 		lq.Factorize(a) | ||||
| 		var l, q Dense | ||||
| 		q.QFromLQ(lq) | ||||
| 		lq.QTo(&q) | ||||
|  | ||||
| 		if !isOrthonormal(&q, 1e-10) { | ||||
| 			t.Errorf("Q is not orthonormal: m = %v, n = %v", m, n) | ||||
| 		} | ||||
|  | ||||
| 		l.LFromLQ(lq) | ||||
| 		lq.LTo(&l) | ||||
|  | ||||
| 		var got Dense | ||||
| 		got.Mul(&l, &q) | ||||
|   | ||||
| @@ -51,7 +51,7 @@ func (lu *LU) updateCond(norm float64) { | ||||
| // The LU factorization is computed with pivoting, and so really the decomposition | ||||
| // is a PLU decomposition where P is a permutation matrix. The individual matrix | ||||
| // factors can be extracted from the factorization using the Permutation method | ||||
| // on Dense, and the LFrom and UFrom methods on TriDense. | ||||
| // on Dense, and the LU LTo and UTo methods. | ||||
| func (lu *LU) Factorize(a Matrix) { | ||||
| 	r, c := a.Dims() | ||||
| 	if r != c { | ||||
| @@ -204,30 +204,30 @@ func (lu *LU) RankOne(orig *LU, alpha float64, x, y *Vector) { | ||||
| 	lu.updateCond(-1) | ||||
| } | ||||
|  | ||||
| // LFromLU extracts the lower triangular matrix from an LU factorization. | ||||
| func (t *TriDense) LFromLU(lu *LU) { | ||||
| // LTo extracts the lower triangular matrix from an LU factorization. | ||||
| func (lu *LU) LTo(dst *TriDense) { | ||||
| 	_, n := lu.lu.Dims() | ||||
| 	t.reuseAs(n, false) | ||||
| 	dst.reuseAs(n, false) | ||||
| 	// Extract the lower triangular elements. | ||||
| 	for i := 0; i < n; i++ { | ||||
| 		for j := 0; j < i; j++ { | ||||
| 			t.mat.Data[i*t.mat.Stride+j] = lu.lu.mat.Data[i*lu.lu.mat.Stride+j] | ||||
| 			dst.mat.Data[i*dst.mat.Stride+j] = lu.lu.mat.Data[i*lu.lu.mat.Stride+j] | ||||
| 		} | ||||
| 	} | ||||
| 	// Set ones on the diagonal. | ||||
| 	for i := 0; i < n; i++ { | ||||
| 		t.mat.Data[i*t.mat.Stride+i] = 1 | ||||
| 		dst.mat.Data[i*dst.mat.Stride+i] = 1 | ||||
| 	} | ||||
| } | ||||
|  | ||||
| // UFromLU extracts the upper triangular matrix from an LU factorization. | ||||
| func (t *TriDense) UFromLU(lu *LU) { | ||||
| // UTo extracts the upper triangular matrix from an LU factorization. | ||||
| func (lu *LU) UTo(dst *TriDense) { | ||||
| 	_, n := lu.lu.Dims() | ||||
| 	t.reuseAs(n, true) | ||||
| 	dst.reuseAs(n, true) | ||||
| 	// Extract the upper triangular elements. | ||||
| 	for i := 0; i < n; i++ { | ||||
| 		for j := i; j < n; j++ { | ||||
| 			t.mat.Data[i*t.mat.Stride+j] = lu.lu.mat.Data[i*lu.lu.mat.Stride+j] | ||||
| 			dst.mat.Data[i*dst.mat.Stride+j] = lu.lu.mat.Data[i*lu.lu.mat.Stride+j] | ||||
| 		} | ||||
| 	} | ||||
| } | ||||
|   | ||||
| @@ -20,12 +20,12 @@ func TestLUD(t *testing.T) { | ||||
| 		var want Dense | ||||
| 		want.Clone(a) | ||||
|  | ||||
| 		lu := &LU{} | ||||
| 		var lu LU | ||||
| 		lu.Factorize(a) | ||||
|  | ||||
| 		var l, u TriDense | ||||
| 		l.LFromLU(lu) | ||||
| 		u.UFromLU(lu) | ||||
| 		lu.LTo(&l) | ||||
| 		lu.UTo(&u) | ||||
| 		var p Dense | ||||
| 		pivot := lu.Pivot(nil) | ||||
| 		p.Permutation(n, pivot) | ||||
| @@ -93,8 +93,8 @@ func TestLURankOne(t *testing.T) { | ||||
| // luReconstruct reconstructs the original A matrix from an LU decomposition. | ||||
| func luReconstruct(lu *LU) *Dense { | ||||
| 	var L, U TriDense | ||||
| 	L.LFromLU(lu) | ||||
| 	U.UFromLU(lu) | ||||
| 	lu.LTo(&L) | ||||
| 	lu.UTo(&U) | ||||
| 	var P Dense | ||||
| 	pivot := lu.Pivot(nil) | ||||
| 	P.Permutation(len(pivot), pivot) | ||||
|   | ||||
| @@ -37,7 +37,7 @@ func (qr *QR) updateCond() { | ||||
| // | ||||
| // The QR decomposition is a factorization of the matrix A such that A = Q * R. | ||||
| // The matrix Q is an orthonormal m×m matrix, and R is an m×n upper triangular matrix. | ||||
| // Q and R can be extracted from the QFromQR and RFromQR methods on Dense. | ||||
| // Q and R can be extracted using the QTo and RTo methods. | ||||
| func (qr *QR) Factorize(a Matrix) { | ||||
| 	m, n := a.Dims() | ||||
| 	if m < n { | ||||
| @@ -60,10 +60,10 @@ func (qr *QR) Factorize(a Matrix) { | ||||
| // TODO(btracey): Add in the "Reduced" forms for extracting the n×n orthogonal | ||||
| // and upper triangular matrices. | ||||
|  | ||||
| // RFromQR extracts the m×n upper trapezoidal matrix from a QR decomposition. | ||||
| func (m *Dense) RFromQR(qr *QR) { | ||||
| // RTo extracts the m×n upper trapezoidal matrix from a QR decomposition. | ||||
| func (qr *QR) RTo(dst *Dense) { | ||||
| 	r, c := qr.qr.Dims() | ||||
| 	m.reuseAs(r, c) | ||||
| 	dst.reuseAs(r, c) | ||||
|  | ||||
| 	// Disguise the QR as an upper triangular | ||||
| 	t := &TriDense{ | ||||
| @@ -76,29 +76,29 @@ func (m *Dense) RFromQR(qr *QR) { | ||||
| 		}, | ||||
| 		cap: qr.qr.capCols, | ||||
| 	} | ||||
| 	m.Copy(t) | ||||
| 	dst.Copy(t) | ||||
|  | ||||
| 	// Zero below the triangular. | ||||
| 	for i := r; i < c; i++ { | ||||
| 		zero(m.mat.Data[i*m.mat.Stride : i*m.mat.Stride+c]) | ||||
| 		zero(dst.mat.Data[i*dst.mat.Stride : i*dst.mat.Stride+c]) | ||||
| 	} | ||||
| } | ||||
|  | ||||
| // QFromQR extracts the m×m orthonormal matrix Q from a QR decomposition. | ||||
| func (m *Dense) QFromQR(qr *QR) { | ||||
| // QTo extracts the m×m orthonormal matrix Q from a QR decomposition. | ||||
| func (qr *QR) QTo(dst *Dense) { | ||||
| 	r, _ := qr.qr.Dims() | ||||
| 	m.reuseAsZeroed(r, r) | ||||
| 	dst.reuseAsZeroed(r, r) | ||||
|  | ||||
| 	// Set Q = I. | ||||
| 	for i := 0; i < r*r; i += r + 1 { | ||||
| 		m.mat.Data[i] = 1 | ||||
| 		dst.mat.Data[i] = 1 | ||||
| 	} | ||||
|  | ||||
| 	// Construct Q from the elementary reflectors. | ||||
| 	work := make([]float64, 1) | ||||
| 	lapack64.Ormqr(blas.Left, blas.NoTrans, qr.qr.mat, qr.tau, m.mat, work, -1) | ||||
| 	lapack64.Ormqr(blas.Left, blas.NoTrans, qr.qr.mat, qr.tau, dst.mat, work, -1) | ||||
| 	work = make([]float64, int(work[0])) | ||||
| 	lapack64.Ormqr(blas.Left, blas.NoTrans, qr.qr.mat, qr.tau, m.mat, work, len(work)) | ||||
| 	lapack64.Ormqr(blas.Left, blas.NoTrans, qr.qr.mat, qr.tau, dst.mat, work, len(work)) | ||||
| } | ||||
|  | ||||
| // SolveQR finds a minimum-norm solution to a system of linear equations defined | ||||
|   | ||||
| @@ -30,16 +30,16 @@ func TestQR(t *testing.T) { | ||||
| 		var want Dense | ||||
| 		want.Clone(a) | ||||
|  | ||||
| 		qr := &QR{} | ||||
| 		var qr QR | ||||
| 		qr.Factorize(a) | ||||
| 		var q, r Dense | ||||
| 		q.QFromQR(qr) | ||||
| 		qr.QTo(&q) | ||||
|  | ||||
| 		if !isOrthonormal(&q, 1e-10) { | ||||
| 			t.Errorf("Q is not orthonormal: m = %v, n = %v", m, n) | ||||
| 		} | ||||
|  | ||||
| 		r.RFromQR(qr) | ||||
| 		qr.RTo(&r) | ||||
|  | ||||
| 		var got Dense | ||||
| 		got.Mul(&q, &r) | ||||
|   | ||||
| @@ -138,42 +138,42 @@ func (svd *SVD) Values(s []float64) []float64 { | ||||
| 	return s | ||||
| } | ||||
|  | ||||
| // UFromSVD extracts the matrix U from the singular value decomposition, storing | ||||
| // the result in-place into the receiver. U is size m×m if svd.Kind() == SVDFull, | ||||
| // of size m×min(m,n) if svd.Kind() == SVDThin, and UFromSVD panics otherwise. | ||||
| func (m *Dense) UFromSVD(svd *SVD) { | ||||
| // UTo extracts the matrix U from the singular value decomposition, storing | ||||
| // the result in-place into dst. U is size m×m if svd.Kind() == SVDFull, | ||||
| // of size m×min(m,n) if svd.Kind() == SVDThin, and UTo panics otherwise. | ||||
| func (svd *SVD) UTo(dst *Dense) { | ||||
| 	kind := svd.kind | ||||
| 	if kind != matrix.SVDFull && kind != matrix.SVDThin { | ||||
| 		panic("mat64: improper SVD kind") | ||||
| 	} | ||||
| 	r := svd.u.Rows | ||||
| 	c := svd.u.Cols | ||||
| 	m.reuseAs(r, c) | ||||
| 	dst.reuseAs(r, c) | ||||
|  | ||||
| 	tmp := &Dense{ | ||||
| 		mat:     svd.u, | ||||
| 		capRows: r, | ||||
| 		capCols: c, | ||||
| 	} | ||||
| 	m.Copy(tmp) | ||||
| 	dst.Copy(tmp) | ||||
| } | ||||
|  | ||||
| // VFromSVD extracts the matrix V from the singular value decomposition, storing | ||||
| // the result in-place into the receiver. V is size n×n if svd.Kind() == SVDFull, | ||||
| // of size n×min(m,n) if svd.Kind() == SVDThin, and VFromSVD panics otherwise. | ||||
| func (m *Dense) VFromSVD(svd *SVD) { | ||||
| // VTo extracts the matrix V from the singular value decomposition, storing | ||||
| // the result in-place into dst. V is size n×n if svd.Kind() == SVDFull, | ||||
| // of size n×min(m,n) if svd.Kind() == SVDThin, and VTo panics otherwise. | ||||
| func (svd *SVD) VTo(dst *Dense) { | ||||
| 	kind := svd.kind | ||||
| 	if kind != matrix.SVDFull && kind != matrix.SVDThin { | ||||
| 		panic("mat64: improper SVD kind") | ||||
| 	} | ||||
| 	r := svd.vt.Rows | ||||
| 	c := svd.vt.Cols | ||||
| 	m.reuseAs(c, r) | ||||
| 	dst.reuseAs(c, r) | ||||
|  | ||||
| 	tmp := &Dense{ | ||||
| 		mat:     svd.vt, | ||||
| 		capRows: r, | ||||
| 		capCols: c, | ||||
| 	} | ||||
| 	m.Copy(tmp.T()) | ||||
| 	dst.Copy(tmp.T()) | ||||
| } | ||||
|   | ||||
| @@ -170,8 +170,8 @@ func TestSVD(t *testing.T) { | ||||
|  | ||||
| func extractSVD(svd *SVD) (s []float64, u, v *Dense) { | ||||
| 	var um, vm Dense | ||||
| 	um.UFromSVD(svd) | ||||
| 	vm.VFromSVD(svd) | ||||
| 	svd.UTo(&um) | ||||
| 	svd.VTo(&vm) | ||||
| 	s = svd.Values(nil) | ||||
| 	return s, &um, &vm | ||||
| } | ||||
|   | ||||
		Reference in New Issue
	
	Block a user
	 kortschak
					kortschak