Golang util.Repartition2x1to3x1函数代码示例

原作者: [db:作者] 来自: [db:来源] 收藏邀请

本文整理汇总了Golang中github.com/hrautila/gomas/util.Repartition2x1to3x1函数的典型用法代码示例。如果您正苦于以下问题：Golang Repartition2x1to3x1函数的具体用法？Golang Repartition2x1to3x1怎么用？Golang Repartition2x1to3x1使用的例子？那么恭喜您, 这里精选的函数代码示例或许可以为您提供帮助。

在下文中一共展示了Repartition2x1to3x1函数的20个代码示例，这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞，您的评价将有助于我们的系统推荐出更棒的Golang代码示例。

示例1: blkReduceTridiagUpper

func blkReduceTridiagUpper(A, tauq, Y, W *cmat.FloatMatrix, lb int, conf *gomas.Config) {
	var ATL, ABR cmat.FloatMatrix
	var A00, A01, A11, A22 cmat.FloatMatrix
	var YT, YB, Y0, Y1, Y2 cmat.FloatMatrix
	var tqT, tqB, tq0, tauq1, tq2 cmat.FloatMatrix
	var v0 float64

	util.Partition2x2(
		&ATL, nil,
		nil, &ABR, A, 0, 0, util.PBOTTOMRIGHT)
	util.Partition2x1(
		&YT,
		&YB, Y, 0, util.PBOTTOM)
	util.Partition2x1(
		&tqT,
		&tqB, tauq, 1, util.PBOTTOM)

	for m(&ATL)-lb > 0 {
		util.Repartition2x2to3x3(&ATL,
			&A00, &A01, nil,
			nil, &A11, nil,
			nil, nil, &A22, A, lb, util.PTOPLEFT)
		util.Repartition2x1to3x1(&YT,
			&Y0,
			&Y1,
			&Y2, Y, lb, util.PTOP)
		util.Repartition2x1to3x1(&tqT,
			&tq0,
			&tauq1,
			&tq2, tauq, lb, util.PTOP)
		// ------------------------------------------------------
		unblkBuildTridiagUpper(&ATL, &tauq1, &YT, W)

		// set subdiagonal entry to unit
		v0 = A01.Get(-1, 0)
		A01.Set(-1, 0, 1.0)

		// A22 := A22 - A01*Y0.T - Y0*A01.T
		blasd.Update2Sym(&A00, &A01, &Y0, -1.0, 1.0, gomas.UPPER, conf)

		// restore subdiagonal entry
		A01.Set(-1, 0, v0)
		// ------------------------------------------------------
		util.Continue3x3to2x2(
			&ATL, nil,
			nil, &ABR, &A00, &A11, &A22, A, util.PTOPLEFT)
		util.Continue3x1to2x1(
			&YT,
			&YB, &Y0, &Y1, Y, util.PTOP)
		util.Continue3x1to2x1(
			&tqT,
			&tqB, &tq0, &tauq1, tauq, util.PTOP)
	}

	if m(&ATL) > 0 {
		unblkReduceTridiagUpper(&ATL, &tqT, W, true)
	}
}

开发者ID:hrautila，项目名称:gomas，代码行数:58，代码来源:trired.go

示例2: blockedQL

/*
 * Blocked QR decomposition with compact WY transform.
 *
 * Compatible with lapack.DGEQRF.
 */
func blockedQL(A, Tvec, Twork, W *cmat.FloatMatrix, lb int, conf *gomas.Config) {
	var ATL, ATR, ABL, ABR, AL cmat.FloatMatrix
	var A00, A01, A10, A11, A22 cmat.FloatMatrix
	var TT, TB cmat.FloatMatrix
	var t0, tau, t2 cmat.FloatMatrix
	var Wrk, w1 cmat.FloatMatrix

	util.Partition2x2(
		&ATL, &ATR,
		&ABL, &ABR, A, 0, 0, util.PBOTTOMRIGHT)
	util.Partition2x1(
		&TT,
		&TB, Tvec, 0, util.PBOTTOM)

	nb := lb
	for m(&ATL)-nb > 0 && n(&ATL)-nb > 0 {
		util.Repartition2x2to3x3(&ATL,
			&A00, &A01, nil,
			&A10, &A11, nil,
			nil, nil, &A22, A, nb, util.PTOPLEFT)
		util.Repartition2x1to3x1(&TT,
			&t0,
			&tau,
			&t2, Tvec, nb, util.PTOP)

		// current block size
		cb, rb := A11.Size()
		if rb < cb {
			cb = rb
		}
		// --------------------------------------------------------
		// decompose righ side AL == /A01\
		//                           \A11/
		w1.SubMatrix(W, 0, 0, cb, 1)
		util.Merge2x1(&AL, &A01, &A11)
		unblockedQL(&AL, &tau, &w1)

		// build block reflector
		unblkQLBlockReflector(Twork, &AL, &tau)

		// update A'tail i.e. A10 and A00 with (I - Y*T*Y.T).T * A'tail
		// compute: C - Y*(C.T*Y*T).T
		ar, ac := A10.Size()
		Wrk.SubMatrix(W, 0, 0, ac, ar)
		updateQLLeft(&A10, &A00, &A11, &A01, Twork, &Wrk, true, conf)
		// --------------------------------------------------------
		util.Continue3x3to2x2(
			&ATL, &ATR,
			&ABL, &ABR, &A00, &A11, &A22, A, util.PTOPLEFT)
		util.Continue3x1to2x1(
			&TT,
			&TB, &t0, &tau, Tvec, util.PTOP)
	}

	// last block with unblocked
	if m(&ATL) > 0 && n(&ATL) > 0 {
		w1.SubMatrix(W, 0, 0, n(&ATL), 1)
		unblockedQL(&ATL, &t0, &w1)
	}
}

开发者ID:hrautila，项目名称:gomas，代码行数:65，代码来源:ql.go

示例3: blockedLQ

/*
 * Blocked LQ decomposition with compact WY transform. As implemented
 * in lapack.DGELQF subroutine.
 */
func blockedLQ(A, Tvec, Twork, W *cmat.FloatMatrix, lb int, conf *gomas.Config) {
	var ATL, ATR, ABL, ABR, AR cmat.FloatMatrix
	var A00, A11, A12, A21, A22 cmat.FloatMatrix
	var TT, TB cmat.FloatMatrix
	var t0, tau, t2 cmat.FloatMatrix
	var Wrk, w1 cmat.FloatMatrix

	util.Partition2x2(
		&ATL, &ATR,
		&ABL, &ABR, A, 0, 0, util.PTOPLEFT)
	util.Partition2x1(
		&TT,
		&TB, Tvec, 0, util.PTOP)

	//nb := conf.LB
	for m(&ABR)-lb > 0 && n(&ABR)-lb > 0 {
		util.Repartition2x2to3x3(&ATL,
			&A00, nil, nil,
			nil, &A11, &A12,
			nil, &A21, &A22, A, lb, util.PBOTTOMRIGHT)
		util.Repartition2x1to3x1(&TT,
			&t0,
			&tau,
			&t2, Tvec, lb, util.PBOTTOM)

		// current block size
		cb, rb := A11.Size()
		if rb < cb {
			cb = rb
		}
		// --------------------------------------------------------
		// decompose left side AL == /A11\
		//                           \A21/
		w1.SubMatrix(W, 0, 0, cb, 1)
		util.Merge1x2(&AR, &A11, &A12)
		unblockedLQ(&AR, &tau, &w1)

		// build block reflector
		unblkBlockReflectorLQ(Twork, &AR, &tau)

		// update A'tail i.e. A21 and A22 with A'*(I - Y*T*Y.T).T
		// compute: C - Y*(C.T*Y*T).T
		ar, ac := A21.Size()
		Wrk.SubMatrix(W, 0, 0, ar, ac)
		updateRightLQ(&A21, &A22, &A11, &A12, Twork, &Wrk, true, conf)
		// --------------------------------------------------------
		util.Continue3x3to2x2(
			&ATL, &ATR,
			&ABL, &ABR, &A00, &A11, &A22, A, util.PBOTTOMRIGHT)
		util.Continue3x1to2x1(
			&TT,
			&TB, &t0, &tau, Tvec, util.PBOTTOM)
	}

	// last block with unblocked
	if m(&ABR) > 0 && n(&ABR) > 0 {
		w1.SubMatrix(W, 0, 0, m(&ABR), 1)
		unblockedLQ(&ABR, &t2, &w1)
	}
}

开发者ID:hrautila，项目名称:gomas，代码行数:64，代码来源:lq.go

示例4: unblkReduceTridiagUpper

/*
 * Reduce upper triangular matrix to tridiagonal.
 *
 * Elementary reflectors Q = H(n-1)...H(2)H(1) are stored on upper
 * triangular part of A. Reflector H(n-1) saved at column A(n) and
 * scalar multiplier to tau[n-1]. If parameter `tail` is true then
 * this function is used to reduce tail part of partially reduced
 * matrix and tau-vector partitioning is starting from last position.
 */
func unblkReduceTridiagUpper(A, tauq, W *cmat.FloatMatrix, tail bool) {
	var ATL, ABR cmat.FloatMatrix
	var A00, a01, a11, A22 cmat.FloatMatrix
	var tqT, tqB, tq0, tauq1, tq2 cmat.FloatMatrix
	var y21 cmat.FloatMatrix
	var v0 float64

	toff := 1
	if tail {
		toff = 0
	}
	util.Partition2x2(
		&ATL, nil,
		nil, &ABR, A, 0, 0, util.PBOTTOMRIGHT)
	util.Partition2x1(
		&tqT,
		&tqB, tauq, toff, util.PBOTTOM)

	for n(&ATL) > 0 {
		util.Repartition2x2to3x3(&ATL,
			&A00, &a01, nil,
			nil, &a11, nil,
			nil, nil, &A22, A, 1, util.PTOPLEFT)
		util.Repartition2x1to3x1(&tqT,
			&tq0,
			&tauq1,
			&tq2, tauq, 1, util.PTOP)
		// set temp vectors for this round
		y21.SetBuf(n(&A00), 1, n(&A00), W.Data())
		// ------------------------------------------------------

		// Compute householder to zero super-diagonal entries
		computeHouseholderRev(&a01, &tauq1)
		tauqv := tauq1.Get(0, 0)

		// set superdiagonal to unit
		v0 = a01.Get(-1, 0)
		a01.Set(-1, 0, 1.0)

		// y21 := A22*a12t
		blasd.MVMultSym(&y21, &A00, &a01, tauqv, 0.0, gomas.UPPER)
		// beta := tauq*a12t*y21
		beta := tauqv * blasd.Dot(&a01, &y21)
		// y21  := y21 - 0.5*beta*a125
		blasd.Axpy(&y21, &a01, -0.5*beta)
		// A22 := A22 - a12t*y21.T - y21*a12.T
		blasd.MVUpdate2Sym(&A00, &a01, &y21, -1.0, gomas.UPPER)

		// restore superdiagonal value
		a01.Set(-1, 0, v0)
		// ------------------------------------------------------
		util.Continue3x3to2x2(
			&ATL, nil,
			nil, &ABR, &A00, &a11, &A22, A, util.PTOPLEFT)
		util.Continue3x1to2x1(
			&tqT,
			&tqB, &tq0, &tauq1, tauq, util.PTOP)
	}
}

开发者ID:hrautila，项目名称:gomas，代码行数:68，代码来源:trired.go

示例5: unblkBlockReflectorLQ

/*
 * like LAPACK/dlafrt.f
 *
 * Build block reflector T from HH reflector stored in TriLU(A) and coefficients
 * in tau.
 *
 * Q = I - Y*T*Y.T; Householder H = I - tau*v*v.T
 *
 * T = | T  z |   z = -tau*T*Y.T*v
 *     | 0  c |   c = tau
 *
 * Q = H(1)H(2)...H(k) building forward here.
 */
func unblkBlockReflectorLQ(T, A, tau *cmat.FloatMatrix) {
	var ATL, ATR, ABL, ABR cmat.FloatMatrix
	var A00, a01, A02, a11, a12, A22 cmat.FloatMatrix
	var TTL, TTR, TBL, TBR cmat.FloatMatrix
	var T00, t01, T02, t11, t12, T22 cmat.FloatMatrix
	var tT, tB cmat.FloatMatrix
	var t0, tau1, t2 cmat.FloatMatrix

	util.Partition2x2(
		&ATL, &ATR,
		&ABL, &ABR, A, 0, 0, util.PTOPLEFT)
	util.Partition2x2(
		&TTL, &TTR,
		&TBL, &TBR, T, 0, 0, util.PTOPLEFT)
	util.Partition2x1(
		&tT,
		&tB, tau, 0, util.PTOP)

	for m(&ABR) > 0 && n(&ABR) > 0 {
		util.Repartition2x2to3x3(&ATL,
			&A00, &a01, &A02,
			nil, &a11, &a12,
			nil, nil, &A22, A, 1, util.PBOTTOMRIGHT)
		util.Repartition2x2to3x3(&TTL,
			&T00, &t01, &T02,
			nil, &t11, &t12,
			nil, nil, &T22, T, 1, util.PBOTTOMRIGHT)
		util.Repartition2x1to3x1(&tT,
			&t0,
			&tau1,
			&t2, tau, 1, util.PBOTTOM)
		// --------------------------------------------------

		// t11 := tau
		tauval := tau1.Get(0, 0)
		if tauval != 0.0 {
			t11.Set(0, 0, tauval)

			// t01 := -tauval*(a01 + A02*a12)
			blasd.Axpby(&t01, &a01, 1.0, 0.0)
			blasd.MVMult(&t01, &A02, &a12, -tauval, -tauval, gomas.NONE)
			// t01 := T00*t01
			blasd.MVMultTrm(&t01, &T00, 1.0, gomas.UPPER)
		}

		// --------------------------------------------------
		util.Continue3x3to2x2(
			&ATL, &ATR,
			&ABL, &ABR, &A00, &a11, &A22, A, util.PBOTTOMRIGHT)
		util.Continue3x3to2x2(
			&TTL, &TTR,
			&TBL, &TBR, &T00, &t11, &T22, T, util.PBOTTOMRIGHT)
		util.Continue3x1to2x1(
			&tT,
			&tB, &t0, &tau1, tau, util.PBOTTOM)
	}
}

开发者ID:hrautila，项目名称:gomas，代码行数:70，代码来源:lq.go

示例6: unblkBlockReflectorRQ

/*
 * like LAPACK/dlafrt.f
 *
 * Build block reflector T from HH reflector stored in TriLU(A) and coefficients
 * in tau.
 *
 * Q = I - Y*T*Y.T; Householder H = I - tau*v*v.T
 *
 * T = | T  0 |   z = -tau*T*Y.T*v
 *     | z  c |   c = tau
 *
 * Q = H(1)H(2)...H(k) building forward here.
 */
func unblkBlockReflectorRQ(T, A, tau *cmat.FloatMatrix) {
	var ATL, ABR cmat.FloatMatrix
	var A00, a10, A20, a11, a21, A22 cmat.FloatMatrix
	var TTL, TBR cmat.FloatMatrix
	var T00, t11, t21, T22 cmat.FloatMatrix
	var tT, tB cmat.FloatMatrix
	var t0, tau1, t2 cmat.FloatMatrix

	util.Partition2x2(
		&ATL, nil,
		nil, &ABR /**/, A, 0, 0, util.PBOTTOMRIGHT)
	util.Partition2x2(
		&TTL, nil,
		nil, &TBR /**/, T, 0, 0, util.PBOTTOMRIGHT)
	util.Partition2x1(
		&tT,
		&tB /**/, tau, 0, util.PBOTTOM)

	for m(&ATL) > 0 && n(&ATL) > 0 {
		util.Repartition2x2to3x3(&ATL,
			&A00, nil, nil,
			&a10, &a11, nil,
			&A20, &a21, &A22 /**/, A, 1, util.PTOPLEFT)
		util.Repartition2x2to3x3(&TTL,
			&T00, nil, nil,
			nil, &t11, nil,
			nil, &t21, &T22 /**/, T, 1, util.PTOPLEFT)
		util.Repartition2x1to3x1(&tT,
			&t0,
			&tau1,
			&t2 /**/, tau, 1, util.PTOP)
		// --------------------------------------------------

		// t11 := tau
		tauval := tau1.Get(0, 0)
		if tauval != 0.0 {
			t11.Set(0, 0, tauval)

			// t21 := -tauval*(a21 + A20*a10)
			blasd.Axpby(&t21, &a21, 1.0, 0.0)
			blasd.MVMult(&t21, &A20, &a10, -tauval, -tauval, gomas.NONE)
			// t21 := T22*t21
			blasd.MVMultTrm(&t21, &T22, 1.0, gomas.LOWER)
		}

		// --------------------------------------------------
		util.Continue3x3to2x2(
			&ATL, nil,
			nil, &ABR /**/, &A00, &a11, &A22, A, util.PTOPLEFT)
		util.Continue3x3to2x2(
			&TTL, nil,
			nil, &TBR /**/, &T00, &t11, &T22, T, util.PTOPLEFT)
		util.Continue3x1to2x1(
			&tT,
			&tB /**/, &t0, &tau1, tau, util.PTOP)
	}
}

开发者ID:hrautila，项目名称:gomas，代码行数:70，代码来源:rq.go

示例7: unblkQLBlockReflector

/*
 * like LAPACK/dlafrt.f
 *
 * Build block reflector T from HH reflector stored in TriLU(A) and coefficients
 * in tau.
 *
 * Q = I - Y*T*Y.T; Householder H = I - tau*v*v.T
 *
 * T = | T  z |   z = -tau*T*Y.T*v
 *     | 0  c |   c = tau
 *
 * Q = H(1)H(2)...H(k) building forward here.
 */
func unblkQLBlockReflector(T, A, tau *cmat.FloatMatrix) {
	var ATL, ABR cmat.FloatMatrix
	var A00, a01, a11, A02, a12, A22 cmat.FloatMatrix
	var TTL, TBR cmat.FloatMatrix
	var T00, t11, t21, T22 cmat.FloatMatrix
	var tT, tB cmat.FloatMatrix
	var t0, tau1, t2 cmat.FloatMatrix

	util.Partition2x2(
		&ATL, nil,
		nil, &ABR, A, 0, 0, util.PBOTTOMRIGHT)
	util.Partition2x2(
		&TTL, nil,
		nil, &TBR, T, 0, 0, util.PBOTTOMRIGHT)
	util.Partition2x1(
		&tT,
		&tB, tau, 0, util.PBOTTOM)

	for m(&ATL) > 0 && n(&ATL) > 0 {
		util.Repartition2x2to3x3(&ATL,
			&A00, &a01, &A02,
			nil, &a11, &a12,
			nil, nil, &A22, A, 1, util.PTOPLEFT)
		util.Repartition2x2to3x3(&TTL,
			&T00, nil, nil,
			nil, &t11, nil,
			nil, &t21, &T22, T, 1, util.PTOPLEFT)
		util.Repartition2x1to3x1(&tT,
			&t0,
			&tau1,
			&t2, tau, 1, util.PTOP)
		// --------------------------------------------------

		// t11 := tau
		tauval := tau1.Get(0, 0)
		if tauval != 0.0 {
			t11.Set(0, 0, tauval)

			// t21 := -tauval*(a12.T + &A02.T*a12)
			blasd.Axpby(&t21, &a12, 1.0, 0.0)
			blasd.MVMult(&t21, &A02, &a01, -tauval, -tauval, gomas.TRANSA)
			// t21 := T22*t01
			blasd.MVMultTrm(&t21, &T22, 1.0, gomas.LOWER)
		}

		// --------------------------------------------------
		util.Continue3x3to2x2(
			&ATL, nil,
			nil, &ABR, &A00, &a11, &A22, A, util.PTOPLEFT)
		util.Continue3x3to2x2(
			&TTL, nil,
			nil, &TBR, &T00, &t11, &T22, T, util.PTOPLEFT)
		util.Continue3x1to2x1(
			&tT,
			&tB, &t0, &tau1, tau, util.PTOP)
	}
}

开发者ID:hrautila，项目名称:gomas，代码行数:70，代码来源:ql.go

示例8: blockedRQ

/*
 * Blocked RQ decomposition with compact WY transform. As implemented
 * in lapack.DGERQF subroutine.
 */
func blockedRQ(A, Tvec, Twork, W *cmat.FloatMatrix, lb int, conf *gomas.Config) {
	var ATL, ABR, AL cmat.FloatMatrix
	var A00, A01, A10, A11, A22 cmat.FloatMatrix
	var TT, TB cmat.FloatMatrix
	var t0, tau, t2 cmat.FloatMatrix
	var Wrk, w1 cmat.FloatMatrix

	util.Partition2x2(
		&ATL, nil,
		nil, &ABR /**/, A, 0, 0, util.PBOTTOMRIGHT)
	util.Partition2x1(
		&TT,
		&TB /**/, Tvec, 0, util.PBOTTOM)

	for m(&ATL)-lb > 0 && n(&ATL)-lb > 0 {
		util.Repartition2x2to3x3(&ATL,
			&A00, &A01, nil,
			&A10, &A11, nil,
			nil, nil, &A22 /**/, A, lb, util.PTOPLEFT)
		util.Repartition2x1to3x1(&TT,
			&t0,
			&tau,
			&t2 /**/, Tvec, n(&A11), util.PTOP)

		// current block size
		cb, rb := A11.Size()
		if rb < cb {
			cb = rb
		}
		// --------------------------------------------------------
		// decompose left side AL == ( A10 A11 )
		w1.SubMatrix(W, 0, 0, cb, 1)
		util.Merge1x2(&AL, &A10, &A11)
		unblockedRQ(&AL, &tau, &w1)

		// build block reflector
		unblkBlockReflectorRQ(Twork, &AL, &tau)

		// compute: (A00 A01)(I - Y*T*Y.T)
		ar, ac := A01.Size()
		Wrk.SubMatrix(W, 0, 0, ar, ac)
		updateRightRQ(&A01, &A00, &A11, &A10, Twork, &Wrk, false, conf)
		// --------------------------------------------------------
		util.Continue3x3to2x2(
			&ATL, nil,
			nil, &ABR, &A00, &A11, &A22, A, util.PTOPLEFT)
		util.Continue3x1to2x1(
			&TT,
			&TB, &t0, &tau, Tvec, util.PTOP)
	}

	// last block with unblocked
	if m(&ATL) > 0 && n(&ATL) > 0 {
		w1.SubMatrix(W, 0, 0, m(&ATL), 1)
		unblockedRQ(&ATL, &TT, &w1)
	}
}

开发者ID:hrautila，项目名称:gomas，代码行数:61，代码来源:rq.go

示例9: unblkReduceTridiagLower

/*
 * Tridiagonal reduction of LOWER triangular symmetric matrix, zero elements below 1st
 * subdiagonal:
 *
 *   A =  (1 - tau*u*u.t)*A*(1 - tau*u*u.T)
 *     =  (I - tau*( 0   0   )) (a11 a12) (I - tau*( 0  0   ))
 *        (        ( 0  u*u.t)) (a21 A22) (        ( 0 u*u.t))
 *
 *  a11, a12, a21 not affected
 *
 *  from LEFT:
 *    A22 = A22 - tau*u*u.T*A22
 *  from RIGHT:
 *    A22 = A22 - tau*A22*u.u.T
 *
 *  LEFT and RIGHT:
 *    A22   = A22 - tau*u*u.T*A22 - tau*(A22 - tau*u*u.T*A22)*u*u.T
 *          = A22 - tau*u*u.T*A22 - tau*A22*u*u.T + tau*tau*u*u.T*A22*u*u.T
 *    [x    = tau*A22*u (vector)]  (SYMV)
 *    A22   = A22 - u*x.T - x*u.T + tau*u*u.T*x*u.T
 *    [beta = tau*u.T*x (scalar)]  (DOT)
 *          = A22 - u*x.T - x*u.T + beta*u*u.T
 *          = A22 - u*(x - 0.5*beta*u).T - (x - 0.5*beta*u)*u.T
 *    [w    = x - 0.5*beta*u]      (AXPY)
 *          = A22 - u*w.T - w*u.T  (SYR2)
 *
 * Result of reduction for N = 5:
 *    ( d  .  .  . . )
 *    ( e  d  .  . . )
 *    ( v1 e  d  . . )
 *    ( v1 v2 e  d . )
 *    ( v1 v2 v3 e d )
 */
func unblkReduceTridiagLower(A, tauq, W *cmat.FloatMatrix) {
	var ATL, ABR cmat.FloatMatrix
	var A00, a11, a21, A22 cmat.FloatMatrix
	var tqT, tqB, tq0, tauq1, tq2 cmat.FloatMatrix
	var y21 cmat.FloatMatrix
	var v0 float64

	util.Partition2x2(
		&ATL, nil,
		nil, &ABR, A, 0, 0, util.PTOPLEFT)
	util.Partition2x1(
		&tqT,
		&tqB, tauq, 0, util.PTOP)

	for m(&ABR) > 0 && n(&ABR) > 0 {
		util.Repartition2x2to3x3(&ATL,
			&A00, nil, nil,
			nil, &a11, nil,
			nil, &a21, &A22, A, 1, util.PBOTTOMRIGHT)
		util.Repartition2x1to3x1(&tqT,
			&tq0,
			&tauq1,
			&tq2, tauq, 1, util.PBOTTOM)
		// set temp vectors for this round
		y21.SetBuf(n(&A22), 1, n(&A22), W.Data())
		// ------------------------------------------------------

		// Compute householder to zero subdiagonal entries
		computeHouseholderVec(&a21, &tauq1)
		tauqv := tauq1.Get(0, 0)

		// set subdiagonal to unit
		v0 = a21.Get(0, 0)
		a21.Set(0, 0, 1.0)

		// y21 := tauq*A22*a21
		blasd.MVMultSym(&y21, &A22, &a21, tauqv, 0.0, gomas.LOWER)
		// beta := tauq*a21.T*y21
		beta := tauqv * blasd.Dot(&a21, &y21)
		// y21  := y21 - 0.5*beta*a21
		blasd.Axpy(&y21, &a21, -0.5*beta)
		// A22 := A22 - a21*y21.T - y21*a21.T
		blasd.MVUpdate2Sym(&A22, &a21, &y21, -1.0, gomas.LOWER)

		// restore subdiagonal
		a21.Set(0, 0, v0)
		// ------------------------------------------------------
		util.Continue3x3to2x2(
			&ATL, nil,
			nil, &ABR, &A00, &a11, &A22, A, util.PBOTTOMRIGHT)
		util.Continue3x1to2x1(
			&tqT,
			&tqB, &tq0, &tauq1, tauq, util.PBOTTOM)
	}
}

开发者ID:hrautila，项目名称:gomas，代码行数:88，代码来源:trired.go

示例10: unblkBuildQL

/*
 * Unblocked code for generating M by N matrix Q with orthogonal columns which
 * are defined as the last N columns of the product of K first elementary
 * reflectors.
 *
 * Parameter nk is last nk elementary reflectors that are not used in computing
 * the matrix Q. Parameter mk length of the first unused elementary reflectors
 * First nk columns are zeroed and subdiagonal mk-nk is set to unit.
 *
 * Compatible with lapack.DORG2L subroutine.
 */
func unblkBuildQL(A, Tvec, W *cmat.FloatMatrix, mk, nk int, mayClear bool) {
	var ATL, ATR, ABL, ABR cmat.FloatMatrix
	var A00, a01, a10, a11, a21, A22 cmat.FloatMatrix
	var tT, tB cmat.FloatMatrix
	var t0, tau1, t2, w12, D cmat.FloatMatrix

	// (mk, nk) = (rows, columns) of upper left partition
	util.Partition2x2(
		&ATL, &ATR,
		&ABL, &ABR, A, mk, nk, util.PTOPLEFT)
	util.Partition2x1(
		&tT,
		&tB, Tvec, nk, util.PTOP)

	// zero the left side
	if nk > 0 && mayClear {
		blasd.Scale(&ABL, 0.0)
		blasd.Scale(&ATL, 0.0)
		D.Diag(&ATL, nk-mk)
		blasd.Add(&D, 1.0)
	}

	for m(&ABR) > 0 && n(&ABR) > 0 {
		util.Repartition2x2to3x3(&ATL,
			&A00, &a01, nil,
			&a10, &a11, nil,
			nil, &a21, &A22, A, 1, util.PBOTTOMRIGHT)
		util.Repartition2x1to3x1(&tT,
			&t0,
			&tau1,
			&t2, Tvec, 1, util.PBOTTOM)
		// ------------------------------------------------------
		w12.SubMatrix(W, 0, 0, a10.Len(), 1)
		applyHouseholder2x1(&tau1, &a01, &a10, &A00, &w12, gomas.LEFT)

		blasd.Scale(&a01, -tau1.Get(0, 0))
		a11.Set(0, 0, 1.0-tau1.Get(0, 0))

		// zero bottom elements
		blasd.Scale(&a21, 0.0)
		// ------------------------------------------------------
		util.Continue3x3to2x2(
			&ATL, &ATR,
			&ABL, &ABR, &A00, &a11, &A22, A, util.PBOTTOMRIGHT)
		util.Continue3x1to2x1(
			&tT,
			&tB, &t0, &tau1, Tvec, util.PBOTTOM)
	}
}

开发者ID:hrautila，项目名称:gomas，代码行数:60，代码来源:qlbld.go

示例11: unblkBuildLQ

/*
 * Unblocked code for generating M by N matrix Q with orthogonal columns which
 * are defined as the first N columns of the product of K first elementary
 * reflectors.
 *
 * Parameters nk = n(A)-K, mk = m(A)-K define the initial partitioning of
 * matrix A.
 *
 *  Q = H(k)H(k-1)...H(1)  , 0 < k <= M, where H(i) = I - tau*v*v.T
 *
 * Computation is ordered as H(k)*H(k-1)...*H(1)*I ie. from bottom to top.
 *
 * If k < M rows k+1:M are cleared and diagonal entries [k+1:M,k+1:M] are
 * set to unit. Then the matrix Q is generated by right multiplying elements below
 * of i'th elementary reflector H(i).
 *
 * Compatible to lapack.xORG2L subroutine.
 */
func unblkBuildLQ(A, Tvec, W *cmat.FloatMatrix, mk, nk int, mayClear bool) {
	var ATL, ATR, ABL, ABR cmat.FloatMatrix
	var A00, a10, a11, a12, a21, A22 cmat.FloatMatrix
	var tT, tB cmat.FloatMatrix
	var t0, tau1, t2, w12, D cmat.FloatMatrix

	util.Partition2x2(
		&ATL, &ATR,
		&ABL, &ABR, A, mk, nk, util.PBOTTOMRIGHT)
	util.Partition2x1(
		&tT,
		&tB, Tvec, mk, util.PBOTTOM)

	// zero the bottom part
	if mk > 0 && mayClear {
		blasd.Scale(&ABL, 0.0)
		blasd.Scale(&ABR, 0.0)
		D.Diag(&ABR)
		blasd.Add(&D, 1.0)
	}

	for m(&ATL) > 0 && n(&ATL) > 0 {
		util.Repartition2x2to3x3(&ATL,
			&A00, nil, nil,
			&a10, &a11, &a12,
			nil, &a21, &A22, A, 1, util.PTOPLEFT)
		util.Repartition2x1to3x1(&tT,
			&t0,
			&tau1,
			&t2, Tvec, 1, util.PTOP)
		// ------------------------------------------------------

		w12.SubMatrix(W, 0, 0, a21.Len(), 1)
		applyHouseholder2x1(&tau1, &a12, &a21, &A22, &w12, gomas.RIGHT)

		blasd.Scale(&a12, -tau1.Get(0, 0))
		a11.Set(0, 0, 1.0-tau1.Get(0, 0))

		// zero
		blasd.Scale(&a10, 0.0)
		// ------------------------------------------------------
		util.Continue3x3to2x2(
			&ATL, &ATR,
			&ABL, &ABR, &A00, &a11, &A22, A, util.PTOPLEFT)
		util.Continue3x1to2x1(
			&tT,
			&tB, &t0, &tau1, Tvec, util.PTOP)
	}
}

开发者ID:hrautila，项目名称:gomas，代码行数:67，代码来源:lqbld.go

示例12: unblockedRQ

/*
 * Unblocked RQ decomposition. As implemented
 * in lapack.DGERQ2 subroutine.
 */
func unblockedRQ(A, Tvec, W *cmat.FloatMatrix) {
	var ATL, ABR cmat.FloatMatrix
	var A00, a11, a01, a10, A22 cmat.FloatMatrix
	var tT, tB cmat.FloatMatrix
	var t0, tau1, t2, w12 cmat.FloatMatrix

	util.Partition2x2(
		&ATL, nil,
		nil, &ABR, A, 0, 0, util.PBOTTOMRIGHT)
	util.Partition2x1(
		&tT,
		&tB, Tvec, 0, util.PBOTTOM)

	for m(&ATL) > 0 && n(&ATL) > 0 {
		util.Repartition2x2to3x3(&ATL,
			&A00, &a01, nil,
			&a10, &a11, nil,
			nil, nil, &A22, A, 1, util.PTOPLEFT)
		util.Repartition2x1to3x1(&tT,
			&t0,
			&tau1,
			&t2, Tvec, 1, util.PTOP)
		// ------------------------------------------------------
		computeHouseholder(&a11, &a10, &tau1)

		w12.SubMatrix(W, 0, 0, a01.Len(), 1)
		applyHouseholder2x1(&tau1, &a10, &a01, &A00, &w12, gomas.RIGHT)
		// ------------------------------------------------------
		util.Continue3x3to2x2(
			&ATL, nil,
			nil, &ABR, &A00, &a11, &A22, A, util.PTOPLEFT)
		util.Continue3x1to2x1(
			&tT,
			&tB, &t0, &tau1, Tvec, util.PTOP)
	}
}

开发者ID:hrautila，项目名称:gomas，代码行数:40，代码来源:rq.go

示例13: blkHessGQvdG

/*
 * Blocked version of Hessenberg reduction algorithm as presented in (1). This
 * version uses compact-WY transformation.
 *
 * Some notes:
 *
 * Elementary reflectors stored in [A11; A21].T are not on diagonal of A11. Update of
 * a block aligned with A11; A21 is as follow
 *
 * 1. Update from left Q(k)*C:
 *                                         c0   0                            c0
 * (I - Y*T*Y.T).T*C = C - Y*(C.T*Y)*T.T = C1 - Y1 * (C1.T.Y1+C2.T*Y2)*T.T = C1-Y1*W
 *                                         C2   Y2                           C2-Y2*W
 *
 * where W = (C1.T*Y1+C2.T*Y2)*T.T and first row of C is not affected by update
 *
 * 2. Update from right C*Q(k):
 *                                       0
 * C - C*Y*T*Y.T = c0;C1;C2 - c0;C1;C2 * Y1 *T*(0;Y1;Y2) = c0; C1-W*Y1; C2-W*Y2
 *                                       Y2
 * where  W = (C1*Y1 + C2*Y2)*T and first column of C is not affected
 *
 */
func blkHessGQvdG(A, Tvec, W *cmat.FloatMatrix, nb int, conf *gomas.Config) *gomas.Error {
	var ATL, ATR, ABL, ABR cmat.FloatMatrix
	var A00, A11, A12, A21, A22, A2 cmat.FloatMatrix
	var tT, tB, td cmat.FloatMatrix
	var t0, t1, t2, T cmat.FloatMatrix
	var V, VT, VB /*V0, V1, V2,*/, Y1, Y2, W0 cmat.FloatMatrix

	//fmt.Printf("blkHessGQvdG...\n")
	T.SubMatrix(W, 0, 0, conf.LB, conf.LB)
	V.SubMatrix(W, conf.LB, 0, m(A), conf.LB)
	td.Diag(&T)

	util.Partition2x2(
		&ATL, &ATR,
		&ABL, &ABR, A, 0, 0, util.PTOPLEFT)
	util.Partition2x1(
		&tT,
		&tB, Tvec, 0, util.PTOP)

	for m(&ABR) > nb+1 && n(&ABR) > nb {
		util.Repartition2x2to3x3(&ATL,
			&A00, nil, nil,
			nil, &A11, &A12,
			nil, &A21, &A22, A, nb, util.PBOTTOMRIGHT)
		util.Repartition2x1to3x1(&tT,
			&t0,
			&t1,
			&t2, Tvec, nb, util.PBOTTOM)

		util.Partition2x1(
			&VT,
			&VB, &V, m(&ATL), util.PTOP)
		// ------------------------------------------------------

		unblkBuildHessGQvdG(&ABR, &T, &VB, nil)
		blasd.Copy(&t1, &td)

		// m(Y) == m(ABR)-1, n(Y) == n(A11)
		Y1.SubMatrix(&ABR, 1, 0, n(&A11), n(&A11))
		Y2.SubMatrix(&ABR, 1+n(&A11), 0, m(&A21)-1, n(&A11))

		// [A01; A02] == ATR := ATR*(I - Y*T*Y.T)
		updateHessRightWY(&ATR, &Y1, &Y2, &T, &VT, conf)

		// A2 = [A12; A22].T
		util.Merge2x1(&A2, &A12, &A22)

		// A2 := A2 - VB*T*A21.T
		be := A21.Get(0, -1)
		A21.Set(0, -1, 1.0)
		blasd.MultTrm(&VB, &T, 1.0, gomas.UPPER|gomas.RIGHT)
		blasd.Mult(&A2, &VB, &A21, -1.0, 1.0, gomas.TRANSB, conf)
		A21.Set(0, -1, be)

		// A2 := (I - Y*T*Y.T).T * A2
		W0.SubMatrix(&V, 0, 0, n(&A2), n(&Y2))
		updateHessLeftWY(&A2, &Y1, &Y2, &T, &W0, conf)

		// ------------------------------------------------------
		util.Continue3x3to2x2(
			&ATL, &ATR,
			&ABL, &ABR, &A00, &A11, &A22, A, util.PBOTTOMRIGHT)
		util.Continue3x1to2x1(
			&tT,
			&tB, &t0, &t1, Tvec, util.PBOTTOM)
	}

	if m(&ABR) > 1 {
		// do the rest with unblocked
		util.Merge2x1(&A2, &ATR, &ABR)
		W0.SetBuf(m(A), 1, m(A), W.Data())
		unblkHessGQvdG(&A2, &tB, &W0, m(&ATR))
	}
	return nil
}

开发者ID:hrautila，项目名称:gomas，代码行数:98，代码来源:hess.go

示例14: blkBuildLQ

func blkBuildLQ(A, Tvec, Twork, W *cmat.FloatMatrix, K, lb int, conf *gomas.Config) {
	var ATL, ATR, ABL, ABR, AL cmat.FloatMatrix
	var A00, A10, A11, A12, A21, A22 cmat.FloatMatrix
	var tT, tB cmat.FloatMatrix
	var t0, tau, t2, Wrk, D, T cmat.FloatMatrix

	nk := n(A) - K
	mk := m(A) - K
	uk := K % lb
	util.Partition2x2(
		&ATL, &ATR,
		&ABL, &ABR, A, mk+uk, nk+uk, util.PBOTTOMRIGHT)
	util.Partition2x1(
		&tT,
		&tB, Tvec, mk+uk, util.PBOTTOM)

	// zero the bottom part __CHECK HERE: nk? or mk?
	if nk+uk > 0 {
		blasd.Scale(&ABL, 0.0)
		if uk > 0 {
			// number of reflectors is not multiple of blocking factor
			// do the first part with unblocked code.
			unblkBuildLQ(&ABR, &tB, W, m(&ABR)-uk, n(&ABR)-uk, true)
		} else {
			// blocking factor is multiple of K
			blasd.Scale(&ABR, 0.0)
			D.Diag(&ABR)
			blasd.Add(&D, 1.0)
		}
	}

	for m(&ATL) > 0 && n(&ATL) > 0 {
		util.Repartition2x2to3x3(&ATL,
			&A00, nil, nil,
			&A10, &A11, &A12,
			nil, &A21, &A22, A, lb, util.PTOPLEFT)
		util.Repartition2x1to3x1(&tT,
			&t0,
			&tau,
			&t2, Tvec, lb, util.PTOP)
		// ------------------------------------------------------
		util.Merge1x2(&AL, &A11, &A12)

		// build block reflector
		T.SubMatrix(Twork, 0, 0, n(&A11), n(&A11))
		unblkBlockReflectorLQ(&T, &AL, &tau)

		// update A21 and A22 with (I - Y*T*Y.T) from right
		ar, ac := A21.Size()
		Wrk.SubMatrix(W, 0, 0, ar, ac)
		updateRightLQ(&A21, &A22, &A11, &A12, &T, &Wrk, false, conf)

		// update current block
		unblkBuildLQ(&AL, &tau, W, 0, n(&A12), false)

		// zero top rows
		blasd.Scale(&A10, 0.0)
		// ------------------------------------------------------
		util.Continue3x3to2x2(
			&ATL, &ATR,
			&ABL, &ABR, &A00, &A11, &A22, A, util.PTOPLEFT)
		util.Continue3x1to2x1(
			&tT,
			&tB, &t0, &tau, Tvec, util.PTOP)
	}
}

开发者ID:hrautila，项目名称:gomas，代码行数:66，代码来源:lqbld.go

示例15: blockedMultQRight

/*
 * Blocked version for computing C = C*Q and C = C*Q.T from elementary reflectors
 * and scalar coefficients.
 *
 * Elementary reflectors and scalar coefficients are used to build block reflector T.
 * Matrix C is updated by applying block reflector T using compact WY algorithm.
 */
func blockedMultQRight(C, A, tau, W *cmat.FloatMatrix, flags, nb int, conf *gomas.Config) {
	var ATL, ATR, ABL, ABR, AL cmat.FloatMatrix
	var A00, A10, A11, A20, A21, A22 cmat.FloatMatrix
	var CL, CR, C0, C1, C2 cmat.FloatMatrix
	var tT, tB cmat.FloatMatrix
	var t0, tau1, t2 cmat.FloatMatrix
	var W0, Wrk, Tw, Twork cmat.FloatMatrix

	var Aref *cmat.FloatMatrix
	var pAdir, pAstart, pDir, pStart, pCstart, pCdir util.Direction
	var bsz, cb, mb int

	// partitioning start and direction
	if flags&gomas.TRANS != 0 {
		// from bottom-right to top-left to produce transpose sequence (C*Q.T)
		pAstart = util.PBOTTOMRIGHT
		pAdir = util.PTOPLEFT
		pStart = util.PBOTTOM
		pDir = util.PTOP
		pCstart = util.PRIGHT
		pCdir = util.PLEFT
		mb = imax(0, m(A)-n(A))
		cb = n(C) - n(A)
		Aref = &ATL
	} else {
		// from top-left to bottom-right to produce normal sequence (C*Q)
		pAstart = util.PTOPLEFT
		pAdir = util.PBOTTOMRIGHT
		pStart = util.PTOP
		pDir = util.PBOTTOM
		pCstart = util.PLEFT
		pCdir = util.PRIGHT
		mb = 0
		cb = 0
		Aref = &ABR
	}

	// intermediate reflector at start of workspace
	Twork.SetBuf(nb, nb, nb, W.Data())
	W0.SetBuf(m(C), nb, m(C), W.Data()[Twork.Len():])

	util.Partition2x2(
		&ATL, &ATR,
		&ABL, &ABR, A, mb, 0, pAstart)
	util.Partition1x2(
		&CL, &CR, C, cb, pCstart)
	util.Partition2x1(
		&tT,
		&tB, tau, 0, pStart)

	transpose := flags&gomas.TRANS != 0

	for m(Aref) > 0 && n(Aref) > 0 {
		util.Repartition2x2to3x3(&ATL,
			&A00, nil, nil,
			&A10, &A11, nil,
			&A20, &A21, &A22, A, nb, pAdir)
		util.Repartition2x1to3x1(&tT,
			&t0,
			&tau1,
			&t2, tau, nb, pDir)

		bsz = n(&A11) // C1 block size must match A11
		util.Repartition1x2to1x3(&CL,
			&C0, &C1, &C2, C, bsz, pCdir)
		// --------------------------------------------------------
		// clear & build block reflector from current block
		util.Merge2x1(&AL, &A11, &A21)
		Tw.SubMatrix(&Twork, 0, 0, bsz, bsz)
		blasd.Scale(&Tw, 0.0)
		unblkQRBlockReflector(&Tw, &AL, &tau1)

		// compute: C*Q.T == C - C*(Y*T*Y.T).T = C - C*Y*T.T*Y.T
		//          C*Q   == C - C*Y*T*Y.T
		Wrk.SubMatrix(&W0, 0, 0, m(&C1), bsz)
		updateWithQTRight(&C1, &C2, &A11, &A21, &Tw, &Wrk, transpose, conf)
		// --------------------------------------------------------
		util.Continue3x3to2x2(
			&ATL, &ATR,
			&ABL, &ABR, &A00, &A11, &A22, A, pAdir)
		util.Continue1x3to1x2(
			&CL, &CR, &C0, &C1, C, pCdir)
		util.Continue3x1to2x1(
			&tT,
			&tB, &t0, &tau1, tau, pDir)
	}

}

开发者ID:hrautila，项目名称:gomas，代码行数:95，代码来源:qrmult.go

示例16: unblockedMultQLeft

/*
 * Unblocked algorith for computing C = Q.T*C and C = Q*C.
 *
 * Q = H(1)H(2)...H(k) where elementary reflectors H(i) are stored on i'th column
 * below diagonal in A.
 *
 * Progressing A from top-left to bottom-right i.e from smaller column numbers
 * to larger, produces H(k)...H(2)H(1) == Q.T. and C = Q.T*C
 *
 * Progressing from bottom-right to top-left produces H(1)H(2)...H(k) == Q and C = Q*C
 */
func unblockedMultQLeft(C, A, tau, w *cmat.FloatMatrix, flags int) {
	var ATL, ATR, ABL, ABR cmat.FloatMatrix
	var A00, a10, a11, A20, a21, A22 cmat.FloatMatrix
	var CT, CB, C0, c1t, C2 cmat.FloatMatrix
	var tT, tB cmat.FloatMatrix
	var t0, tau1, t2, w1 cmat.FloatMatrix

	var Aref *cmat.FloatMatrix
	var pAdir, pAstart, pDir, pStart util.Direction
	var mb, tb, nb int

	// partitioning start and direction
	if flags&gomas.TRANS != 0 {
		// from top-left to bottom-right to produce transposed sequence (Q.T*C)
		pAstart = util.PTOPLEFT
		pAdir = util.PBOTTOMRIGHT
		pStart = util.PTOP
		pDir = util.PBOTTOM
		mb = 0
		tb = 0
		nb = 0
		Aref = &ABR
	} else {
		// from bottom-right to top-left to produce normal sequence (Q*C)
		pAstart = util.PBOTTOMRIGHT
		pAdir = util.PTOPLEFT
		pStart = util.PBOTTOM
		pDir = util.PTOP
		mb = imax(0, m(A)-n(A))
		nb = imax(0, n(A)-m(A))
		tb = imax(0, tau.Len()-n(A))
		Aref = &ATL
	}

	util.Partition2x2(
		&ATL, &ATR,
		&ABL, &ABR, A, mb, nb, pAstart)
	util.Partition2x1(
		&CT,
		&CB, C, mb, pStart)
	util.Partition2x1(
		&tT,
		&tB, tau, tb, pStart)

	for m(Aref) > 0 && n(Aref) > 0 {
		util.Repartition2x2to3x3(&ATL,
			&A00, nil, nil,
			&a10, &a11, nil,
			&A20, &a21, &A22, A, 1, pAdir)
		util.Repartition2x1to3x1(&CT,
			&C0,
			&c1t,
			&C2, C, 1, pDir)
		util.Repartition2x1to3x1(&tT,
			&t0,
			&tau1,
			&t2, tau, 1, pDir)

		// --------------------------------------------------------

		w1.SubMatrix(w, 0, 0, c1t.Len(), 1)
		applyHouseholder2x1(&tau1, &a21, &c1t, &C2, &w1, gomas.LEFT)

		// --------------------------------------------------------
		util.Continue3x3to2x2(
			&ATL, &ATR,
			&ABL, &ABR, &A00, &a11, &A22, A, pAdir)
		util.Continue3x1to2x1(
			&CT,
			&CB, &C0, &c1t, C, pDir)
		util.Continue3x1to2x1(
			&tT,
			&tB, &t0, &tau1, tau, pDir)
	}
}

开发者ID:hrautila，项目名称:gomas，代码行数:86，代码来源:qrmult.go

示例17: unblockedMultRQRight

/*
 * Unblocked algorith for computing C = C*Q.T and C = C*Q.
 *
 * Q = H(0)H(1)...H(K-1) where elementary reflectors H(i) are stored on i'th row
 * right of diagonal in A.
 *
 * Progressing A from top-left to bottom-right i.e from smaller column numbers
 * to larger, produces C*H(0)H(1)...H(K-1) == C*Q.
 *
 * Progressing from bottom-right to top-left produces C*H(K-1)...H(1)H(0) == C*Q.T.
 */
func unblockedMultRQRight(C, A, tau, w *cmat.FloatMatrix, flags int) {
	var ATL, ABR cmat.FloatMatrix
	var A00, a10, a11, A22 cmat.FloatMatrix
	var CL, CR, C0, c1, C2 cmat.FloatMatrix
	var tT, tB cmat.FloatMatrix
	var t0, tau1, t2, w1 cmat.FloatMatrix

	var Aref *cmat.FloatMatrix
	var pAdir, pAstart, pDir, pStart, pCstart, pCdir util.Direction
	var cb, mb, tb, nb int

	// partitioning start and direction
	if flags&gomas.TRANS != 0 {
		// from bottom-right to top-left to produce transpose sequence (C*Q.T)
		pAstart = util.PBOTTOMRIGHT
		pAdir = util.PTOPLEFT
		pStart = util.PBOTTOM
		pDir = util.PTOP
		pCstart = util.PRIGHT
		pCdir = util.PLEFT
		mb = 0
		nb = 0
		cb = 0
		tb = 0
		Aref = &ATL
	} else {
		// from top-left to bottom-right to produce normal sequence (C*Q)
		pAstart = util.PTOPLEFT
		pAdir = util.PBOTTOMRIGHT
		pStart = util.PTOP
		pDir = util.PBOTTOM
		pCstart = util.PLEFT
		pCdir = util.PRIGHT
		mb = imax(0, m(A)-n(A))
		nb = imax(0, n(A)-m(A))
		cb = imax(0, n(C)-m(A))
		tb = imax(0, tau.Len()-n(A))
		Aref = &ABR
	}

	util.Partition2x2(
		&ATL, nil,
		nil, &ABR /**/, A, mb, nb, pAstart)
	util.Partition1x2(
		&CL, &CR /**/, C, cb, pCstart)
	util.Partition2x1(
		&tT,
		&tB /**/, tau, tb, pStart)

	w1.SubMatrix(w, 0, 0, m(C), 1)

	for m(Aref) > 0 && n(Aref) > 0 {
		util.Repartition2x2to3x3(&ATL,
			&A00, nil, nil,
			&a10, &a11, nil,
			nil, nil, &A22 /**/, A, 1, pAdir)
		util.Repartition1x2to1x3(&CL,
			&C0, &c1, &C2 /**/, C, 1, pCdir)
		util.Repartition2x1to3x1(&tT,
			&t0,
			&tau1,
			&t2 /**/, tau, 1, pDir)
		// --------------------------------------------------------
		applyHouseholder2x1(&tau1, &a10, &c1, &C0, &w1, gomas.RIGHT)
		// --------------------------------------------------------
		util.Continue3x3to2x2(
			&ATL, nil,
			nil, &ABR /**/, &A00, &a11, &A22, A, pAdir)
		util.Continue1x3to1x2(
			&CL, &CR /**/, &C0, &c1, C, pCdir)
		util.Continue3x1to2x1(
			&tT,
			&tB /**/, &t0, &tau1, tau, pDir)
	}
}

开发者ID:hrautila，项目名称:gomas，代码行数:86，代码来源:rqmult.go

示例18: blockedMultQTLeft

/*
 * Blocked version for computing C = Q*C and C = Q.T*C with block reflector.
 *
 * Block reflector T is [T(0), T(1), ... T(k-1)], conf.LB*n(A) matrix where
 * where each T(n), expect T(k-1), is conf.LB*conf.LB. T(k-1) is IB*IB where
 * IB = imin(LB, K%LB)
 */
func blockedMultQTLeft(C, A, T, W *cmat.FloatMatrix, flags int, conf *gomas.Config) *gomas.Error {
	var ATL, ATR, ABL, ABR cmat.FloatMatrix
	var A00, A10, A11, A20, A21, A22 cmat.FloatMatrix
	var CT, CB, C0, C1, C2 cmat.FloatMatrix
	var TL, TR, T00, T01, T02 cmat.FloatMatrix
	var Wrk cmat.FloatMatrix

	var Aref *cmat.FloatMatrix
	var pAdir, pAstart, pCdir, pCstart, pTstart, pTdir util.Direction
	var bsz, mb, tb int

	lb := conf.LB
	if conf.LB == 0 {
		lb = m(T)
	}
	transpose := flags&gomas.TRANS != 0
	nb := lb
	//W0 := cmat.MakeMatrix(n(C), conf.LB, W.Data())

	// partitioning start and direction
	if flags&gomas.TRANS != 0 {
		// from top-left to bottom-right to produce transposed sequence (Q.T*C)
		pAstart = util.PTOPLEFT
		pAdir = util.PBOTTOMRIGHT
		pCstart = util.PTOP
		pCdir = util.PBOTTOM
		pTstart = util.PLEFT
		pTdir = util.PRIGHT
		mb = 0
		tb = nb
		Aref = &ABR
	} else {
		// from bottom-right to top-left to produce normal sequence (Q*C)
		pAstart = util.PBOTTOMRIGHT
		pAdir = util.PTOPLEFT
		pCstart = util.PBOTTOM
		pCdir = util.PTOP
		pTstart = util.PRIGHT
		pTdir = util.PLEFT
		mb = m(A) - n(A)
		Aref = &ATL
		// if N%LB != 0 then the last T is not of size LB and we need
		// adjust first repartitioning accordingly.
		tb = n(A) % lb
	}

	util.Partition2x2(
		&ATL, &ATR,
		&ABL, &ABR, A, mb, 0, pAstart)
	util.Partition1x2(
		&TL, &TR, T, 0, pTstart)
	util.Partition2x1(
		&CT,
		&CB, C, mb, pCstart)

	nb = tb
	for m(Aref) > 0 && n(Aref) > 0 {
		util.Repartition2x2to3x3(&ATL,
			&A00, nil, nil,
			&A10, &A11, nil,
			&A20, &A21, &A22, A, nb, pAdir)
		util.Repartition1x2to1x3(&TL,
			&T00, &T01, &T02, T, nb, pTdir)

		bsz = n(&A11) // must match A11 block size
		util.Repartition2x1to3x1(&CT,
			&C0,
			&C1,
			&C2, C, bsz, pCdir)
		// --------------------------------------------------------
		tr, tc := T01.Size()
		// compute: Q.T*C == C - Y*(C.T*Y*T).T  transpose == true
		//          Q*C   == C - C*Y*T*Y.T      transpose == false
		Wrk.SubMatrix(W, 0, 0, n(&C1), bsz)
		if tr != tc {
			// this happens when n(A) not multiple of LB
			var Tmp cmat.FloatMatrix
			Tmp.SubMatrix(&T01, 0, 0, tc, tc)
			updateWithQTLeft(&C1, &C2, &A11, &A21, &Tmp, &Wrk, transpose, conf)
		} else {
			updateWithQTLeft(&C1, &C2, &A11, &A21, &T01, &Wrk, transpose, conf)
		}

		// ---------------