cgo: add Dgebal

cgo/lapack.go

@@ -227,6 +227,74 @@ func (impl Implementation) Dpotrf(ul blas.Uplo, n int, a []float64, lda int) (ok
 	return lapacke.Dpotrf(ul, n, a, lda)
 }
 
+// Dgebal balances an n×n matrix A. Balancing consists of two stages, permuting
+// and scaling. Both steps are optional and depend on the value of job.
+//
+// Permuting consists of applying a permutation matrix P such that the matrix
+// that results from P^T*A*P takes the upper block triangular form
+//            [ T1  X  Y  ]
+//  P^T A P = [  0  B  Z  ],
+//            [  0  0  T2 ]
+// where T1 and T2 are upper triangular matrices and B contains at least one
+// nonzero off-diagonal element in each row and column. The indices ilo and ihi
+// mark the starting and ending columns of the submatrix B. The eigenvalues of A
+// isolated in the first 0 to ilo-1 and last ihi+1 to n-1 elements on the
+// diagonal can be read off without any roundoff error.
+//
+// Scaling consists of applying a diagonal similarity transformation D such that
+// D^{-1}*B*D has the 1-norm of each row and its corresponding column nearly
+// equal. The output matrix is
+//  [ T1     X*D          Y    ]
+//  [  0  inv(D)*B*D  inv(D)*Z ].
+//  [  0      0           T2   ]
+// Scaling may reduce the 1-norm of the matrix, and improve the accuracy of
+// the computed eigenvalues and/or eigenvectors.
+//
+// job specifies the operations that will be performed on A.
+// If job is lapack.None, Dgebal sets scale[i] = 1 for all i and returns ilo=0, ihi=n-1.
+// If job is lapack.Permute, only permuting will be done.
+// If job is lapack.Scale, only scaling will be done.
+// If job is lapack.PermuteScale, both permuting and scaling will be done.
+//
+// On return, if job is lapack.Permute or lapack.PermuteScale, it will hold that
+//  A[i,j] == 0,   for i > j and j ∈ {0, ..., ilo-1, ihi+1, ..., n-1}.
+// If job is lapack.None or lapack.Scale, or if n == 0, it will hold that
+//  ilo == 0 and ihi == n-1.
+//
+// On return, scale will contain information about the permutations and scaling
+// factors applied to A. If π(j) denotes the index of the column interchanged
+// with column j, and D[j,j] denotes the scaling factor applied to column j,
+// then
+//  scale[j] == π(j),     for j ∈ {0, ..., ilo-1, ihi+1, ..., n-1},
+//           == D[j,j],   for j ∈ {ilo, ..., ihi}.
+// scale must have length equal to n, otherwise Dgebal will panic.
+//
+// Dgebal is an internal routine. It is exported for testing purposes.
+func (impl Implementation) Dgebal(job lapack.Job, n int, a []float64, lda int, scale []float64) (ilo, ihi int) {
+	switch job {
+	default:
+		panic(badJob)
+	case lapack.None, lapack.Permute, lapack.Scale, lapack.PermuteScale:
+	}
+	checkMatrix(n, n, a, lda)
+	if len(scale) != n {
+		panic("lapack: bad length of scale")
+	}
+
+	ilo32 := make([]int32, 1)
+	ihi32 := make([]int32, 1)
+	lapacke.Dgebal(job, n, a, lda, ilo32, ihi32, scale)
+	ilo = int(ilo32[0]) - 1
+	ihi = int(ihi32[0]) - 1
+	for j := 0; j < ilo; j++ {
+		scale[j]--
+	}
+	for j := ihi + 1; j < n; j++ {
+		scale[j]--
+	}
+	return ilo, ihi
+}
+
 // Dbdsqr performs a singular value decomposition of a real n×n bidiagonal matrix.
 //
 // The SVD of the bidiagonal matrix B is

cgo/lapack_test.go

@@ -64,6 +64,10 @@ func TestDpotrf(t *testing.T) {
 	testlapack.DpotrfTest(t, impl)
 }
 
+func TestDgebal(t *testing.T) {
+	testlapack.DgebalTest(t, impl)
+}
+
 func TestDgebd2(t *testing.T) {
 	testlapack.Dgebd2Test(t, blockedTranslate{impl})
 }