Golang FloatMatrix.Set方法代码示例

func TestPartition2D(t *testing.T) {
	var ATL, ATR, ABL, ABR, As cmat.FloatMatrix
	var A00, a01, A02, a10, a11, a12, A20, a21, A22 cmat.FloatMatrix

	csource := cmat.NewFloatConstSource(1.0)
	A := cmat.NewMatrix(6, 6)
	As.SubMatrix(A, 1, 1, 4, 4)
	As.SetFrom(csource)

	Partition2x2(&ATL, &ATR, &ABL, &ABR, &As, 0, 0, PTOPLEFT)
	t.Logf("ATL:\n%v\n", &ATL)

	t.Logf("n(ATL)=%d, n(As)=%d\n", n(&ATL), n(&As))
	k := 0
	for n(&ATL) < n(&As) && k < n(&As) {
		Repartition2x2to3x3(&ATL,
			&A00, &a01, &A02,
			&a10, &a11, &a12,
			&A20, &a21, &A22, &As, 1, PBOTTOMRIGHT)
		t.Logf("n(A00)=%d, n(a01)=%d, n(A02)=%d\n", n(&A00), n(&a01), n(&A02))
		t.Logf("n(a10)=%d, n(a11)=%d, n(a12)=%d\n", n(&a10), n(&a11), n(&a12))
		t.Logf("n(A20)=%d, n(a21)=%d, n(A22)=%d\n", n(&A20), n(&a21), n(&A22))
		//t.Logf("n(a12)=%d [%d], n(a11)=%d\n", n(&a12), a12.Len(), a11.Len())
		a11.Set(0, 0, a11.Get(0, 0)+1.0)
		addConst(&a21, -2.0)

		Continue3x3to2x2(&ATL, &ATR, &ABL, &ABR, &A00, &a11, &A22, &As, PBOTTOMRIGHT)
		t.Logf("n(ATL)=%d, n(As)=%d\n", n(&ATL), n(&As))
		k += 1
	}
	t.Logf("A:\n%v\n", A)
}

/*
 * 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)
	}
}

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)
	}
}

/*
 * 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)
	}
}

/*
 * 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)
	}
}

/*
 * 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)
	}
}

/*
 * 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)
	}
}

/*
 * 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)
	}
}

/*
 * 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)
	}
}

/*
 * Apply diagonal pivot (row and column swapped) to symmetric matrix blocks.
 *
 * LOWER triangular; moving from top-left to bottom-right
 *
 *    -----------------------
 *    | d
 *    | x P1 x  x  x  P2     -- current row/col 'srcix'
 *    | x S2 d  x  x  x
 *    | x S2 x  d  x  x
 *    | x S2 x  x  d  x
 *    | x P2 D2 D2 D2 P3     -- swap with row/col 'dstix'
 *    | x S3 x  x  x  D3 d
 *    | x S3 x  x  x  D3 x d
 *         (AR)
 */
func applyBKPivotSymLower(AR *cmat.FloatMatrix, srcix, dstix int) {
	var s, d cmat.FloatMatrix
	// S2 -- D2
	s.SubMatrix(AR, srcix+1, srcix, dstix-srcix-1, 1)
	d.SubMatrix(AR, dstix, srcix+1, 1, dstix-srcix-1)
	blasd.Swap(&s, &d)
	// S3 -- D3
	s.SubMatrix(AR, dstix+1, srcix, m(AR)-dstix-1, 1)
	d.SubMatrix(AR, dstix+1, dstix, m(AR)-dstix-1, 1)
	blasd.Swap(&s, &d)
	// swap P1 and P3
	p1 := AR.Get(srcix, srcix)
	p3 := AR.Get(dstix, dstix)
	AR.Set(srcix, srcix, p3)
	AR.Set(dstix, dstix, p1)
	return
}

/*
 * Apply diagonal pivot (row and column swapped) to symmetric matrix blocks.
 *
 * UPPER triangular; moving from bottom-right to top-left
 *
 *    d x D3 x  x  x  S3 x |
 *      d D3 x  x  x  S3 x |
 *        P3 D2 D2 D2 P2 x |  -- dstinx
 *           d  x  x  S2 x |
 *              d  x  S2 x |
 *                 d  S2 x |
 *                    P1 x |  -- srcinx
 *                       d |
 *    ----------------------
 *               (ABR)
 */
func applyBKPivotSymUpper(AR *cmat.FloatMatrix, srcix, dstix int) {
	var s, d cmat.FloatMatrix
	// AL is ATL, AR is ATR; P1 is AL[srcix, srcix];
	// S2 -- D2
	s.SubMatrix(AR, dstix+1, srcix, srcix-dstix-1, 1)
	d.SubMatrix(AR, dstix, dstix+1, 1, srcix-dstix-1)
	blasd.Swap(&s, &d)
	// S3 -- D3
	s.SubMatrix(AR, 0, srcix, dstix, 1)
	d.SubMatrix(AR, 0, dstix, dstix, 1)
	blasd.Swap(&s, &d)
	// swap P1 and P3
	p1 := AR.Get(srcix, srcix)
	p3 := AR.Get(dstix, dstix)
	AR.Set(srcix, srcix, p3)
	AR.Set(dstix, dstix, p1)
	return
}

func TestViews(t *testing.T) {
	var As cmat.FloatMatrix
	N := 7
	A := cmat.NewMatrix(N, N)
	zeromean := cmat.NewFloatNormSource()
	A.SetFrom(zeromean)
	dlast := A.Get(-1, -1)
	t.Logf("A[-1,-1] = %10.3e\n", dlast)
	for i := 0; i < N; i++ {
		As.SubMatrix(A, i, i)
		As.Set(0, 0, As.Get(0, 0)+float64(i+1))
		if N < 10 {
			t.Logf("A[%d,%d]   = %10.3e\n", i, i, As.Get(0, 0))
		}
	}
	ok := float64(N)+dlast == A.Get(-1, -1)
	nC := int(A.Get(-1, -1) - dlast)
	t.Logf("add 1.0 %d times [%10.3e to %10.3e]: %v\n", nC, dlast, A.Get(-1, -1), ok)
}

/*
 * Generates the real orthogonal matrix Q which is defined as the product of K elementary
 * reflectors of order N embedded in matrix A as returned by TRDReduce().
 *
 *   A     On entry tridiagonal reduction as returned by TRDReduce().
 *         On exit the orthogonal matrix Q.
 *
 *  tau    Scalar coefficients of elementary reflectors.
 *
 *  W      Workspace
 *
 *  K      Number of reflectors , 0 < K < N
 *
 *  flags  LOWER or UPPER
 *
 *  confs  Optional blocking configuration
 *
 * If flags has UPPER set then
 *    Q = H(K)...H(1)H(0) where 0 < K < N-1
 *
 * If flags has LOWR set then
 *    Q = H(0)H(1)...H(K) where 0 < K < N-1
 */
func TRDBuild(A, tau, W *cmat.FloatMatrix, K, flags int, confs ...*gomas.Config) *gomas.Error {
	var err *gomas.Error = nil
	var Qh, tauh cmat.FloatMatrix
	var s, d cmat.FloatMatrix

	if K > m(A)-1 {
		K = m(A) - 1
	}

	switch flags & (gomas.LOWER | gomas.UPPER) {
	case gomas.LOWER:
		// Shift Q matrix embedded in A right and fill first column
		// unit column vector
		for j := m(A) - 1; j > 0; j-- {
			s.SubMatrix(A, j, j-1, m(A)-j, 1)
			d.SubMatrix(A, j, j, m(A)-j, 1)
			blasd.Copy(&d, &s)
			A.Set(0, j, 0.0)
		}
		// zero first column and set first entry to one
		d.Column(A, 0)
		blasd.Scale(&d, 0.0)
		d.Set(0, 0, 1.0)

		Qh.SubMatrix(A, 1, 1, m(A)-1, m(A)-1)
		tauh.SubMatrix(tau, 0, 0, m(A)-1, 1)
		err = QRBuild(&Qh, &tauh, W, K, confs...)

	case gomas.UPPER:
		// Shift Q matrix embedded in A left and fill last column
		// unit column vector
		for j := 1; j < m(A); j++ {
			s.SubMatrix(A, 0, j, j, 1)
			d.SubMatrix(A, 0, j-1, j, 1)
			blasd.Copy(&d, &s)
			A.Set(-1, j-1, 0.0)
		}
		// zero last column and set last entry to one
		d.Column(A, m(A)-1)
		blasd.Scale(&d, 0.0)
		d.Set(-1, 0, 1.0)

		Qh.SubMatrix(A, 0, 0, m(A)-1, m(A)-1)
		tauh.SubMatrix(tau, 0, 0, m(A)-1, 1)
		err = QLBuild(&Qh, &tauh, W, K, confs...)
	}
	if err != nil {
		err.Update("TRDBuild")
	}
	return err
}

/*
 * Compute A[i:i+nr,c1;c2] = A[i:i+nr,c1;c2]*G(c,s)
 *
 * Applies Givens rotation to nr rows of columns c1, c2 of A starting at row.
 *
 */
func ApplyGivensRight(A *cmat.FloatMatrix, c1, c2, row, nrow int, cos, sin float64) {
	if row >= m(A) {
		return
	}
	if c1 == c2 {
		// one column
		for k := row; k < row+nrow; k++ {
			v0 := A.Get(k, c1)
			A.Set(k, c1, cos*v0)
		}
		return
	}
	for k := row; k < row+nrow; k++ {
		v0 := A.Get(k, c1)
		v1 := A.Get(k, c2)
		y0 := cos*v0 + sin*v1
		y1 := cos*v1 - sin*v0
		A.Set(k, c1, y0)
		A.Set(k, c2, y1)
	}
}

/*
 * Compute A[i:i+1,j:j+nc] = G(c,s)*A[i:i+1,j:j+nc]
 *
 * Applies Givens rotation to nc columns on rows r1,r2 of A starting from col.
 */
func ApplyGivensLeft(A *cmat.FloatMatrix, r1, r2, col, ncol int, cos, sin float64) {
	if col >= n(A) {
		return
	}
	if r1 == r2 {
		// one row
		for k := col; k < col+ncol; k++ {
			v0 := A.Get(r1, k)
			A.Set(r1, k, cos*v0)
		}
		return
	}
	for k := col; k < col+ncol; k++ {
		v0 := A.Get(r1, k)
		v1 := A.Get(r2, k)
		y0 := cos*v0 + sin*v1
		y1 := cos*v1 - sin*v0
		A.Set(r1, k, y0)
		A.Set(r2, k, y1)
	}
}