计算数学
go语言
go语言中，也可以把矩阵和向量定义为结构体。进而可以类似于C++或C#定义其相应的方法。 在数值运算时，矩阵和向量类直接声明为所定义的结构体。不过由于go中的切片较为简单，个人认为没有必要再去单独定义。 下面给出一些简单常用的矩阵和向量函数。
需要注意的是，切片在使用时，是地址引用类型。
package main
import (
"bufio"
"fmt"
"io"
"math"
"os"
"strconv"
"strings"
)
func Vec(N int) (V []float64) {
/*------------------------------------------------------
! Author : Song Yezhi
! verison : 2020-5-24 17:50
!
!
! -----------------------------------------------------
! Input Parameters :
!
! Output Parameters :
!
------------------------------------------------------*/
V = make([]float64, N)
for i := 0; i < N; i++ {
V[i] = 0.0
}
return
}
func Mat(M int, N int) (A [][]float64) {
/*------------------------------------------------------
! Author : Song Yezhi
! verison : 2020-5-23 23:43
!
! creat a two dimention slices (matrix)
! -----------------------------------------------------
! Input Parameters :
!
! Output Parameters :
!
------------------------------------------------------*/
A = make([][]float64, M)
for i := 0; i < M; i++ {
A[i] = make([]float64, N)
}
for i := 0; i < M; i++ {
for j := 0; j < N; j++ {
A[i][j] = 0.0
}
}
return
}
func MatDim(A [][]float64) (M int, N int) {
/*------------------------------------------------------
! Author : Song Yezhi
! verison : 2020-5-24 12:45
------------------------------------------------------*/
M = len(A)
N = len(A[M-1])
return
}
func VecOutput(V []float64) {
N := len(V)
for i := 0; i < N; i++ {
fmt.Printf("%18.6f\n", V[i])
}
}
func VecDot(x1 []float64, x2 []float64) (dot float64) {
/*------------------------------------------------------
! Author : Song Yezhi
! verison : 2020-5-19 22:53
!
! dot product
------------------------------------------------------*/
N := len(x1)
var i int
dot = 0.0
for i = 0; i < N; i++ {
dot = dot + x1[i]*x2[i]
}
return
}
func MatVec(A [][]float64, x []float64) (b []float64) {
/*------------------------------------------------------
! Author : Song Yezhi
! verison : 2020-5-19 22:25
! b = Ax
!
! -----------------------------------------------------
! Input Parameters :
! A[M][N] ----MATRIX
! Output Parameters :
! b[M]-----vector
------------------------------------------------------*/
M := len(A)
N := len(A[0])
//the size of mat
//b = make([]float64, M)
b = Vec(M)
var sum float64
for i := 0; i < M; i++ {
sum = 0.0
for j := 0; j < N; j++ {
sum = sum + A[i][j]*x[j]
}
b[i] = sum
}
return
}
func MatMul(A [][]float64, B [][]float64) (C [][]float64) {
/*------------------------------------------------------
! Author : Song Yezhi
! verison : 2020-5-23 23:27
! C = A*B
!
! -----------------------------------------------------
! Input Parameters :
! A[M][N]
! B[N][P]
! Output Parameters :
! C[M][P]
------------------------------------------------------*/
M := len(A)
N := len(A[M-1])
P := len(B[N-1])
var i, j, k int
C = Mat(M, P)
var sum float64
for i = 0; i < M; i++ {
for j = 0; j < P; j++ {
sum = 0.0
for k = 0; k < N; k++ {
sum = sum + A[i][k]*B[k][j]
C[i][j] = sum
}
}
}
return
}
func MatOutput(A [][]float64) {
/*------------------------------------------------------
! Author : Song Yezhi
! verison : 2020-5-24 16:54
!
------------------------------------------------------*/
M, N := MatDim(A)
for i := 0; i < M; i++ {
for j := 0; j < N; j++ {
fmt.Printf("%18.6f ", A[i][j])
}
fmt.Printf("\n")
}
}
func LDL(A [][]float64) (L [][]float64, D []float64) {
/*------------------------------------------------------
! Author : Song Yezhi
! verison : 2020-5-24 17:24
!
! A=LDL'
! -----------------------------------------------------
! Input Parameters :
! A
! Output Parameters :
! L
! D
------------------------------------------------------*/
M, N := MatDim(A)
if M != N {
fmt.Printf("the rows and cols is not equal")
}
G := Mat(N, N)
L = Mat(N, N)
D = Vec(N)
var i, j, k int
D[0] = A[0][0]
var tmp1 float64
for i = 1; i < N; i++ {
for j = 0; j <= i-1; j++ {
tmp1 = 0.0
for k = 0; k <= j-1; k++ {
tmp1 = tmp1 + G[i][k]*L[j][k]
}
G[i][j] = A[i][j] - tmp1
}
for j = 0; j <= i-1; j++ {
L[i][j] = G[i][j] / D[j]
}
//comput D[i]
tmp1 = 0.0
for k = 0; k <= i-1; k++ {
tmp1 = tmp1 + G[i][k]*L[i][k]
}
D[i] = A[i][i] - tmp1
}
//matOutput(L)
//set the diag elements
for i = 0; i < N; i++ {
L[i][i] = 1.0
}
return
}
func AddNeq(ATA [][]float64, ATb []float64, H []float64, oc float64) {
/*------------------------------------------------------
Created : Song Yezhi 2022-5-21 14:15
add a new observation to normal equations
(ATA+ H'H) x = ATb + H'oc
ATA = ATA + H'H
ATb = ATb + H'oc
--------------------------------------------------------
Input Parameters :
ATA
ATb
H
oc
Output Parameters :
ATA --- updated
ATb --- updated
--------------------------------------------------------
Email : song.yz@foxmail.com
Copyrigt (C) : Chinese Academy of Sciences
All rights reserved, 2022
-------------------------------------------------------*/
N := len(ATA)
var i, j int
for i = 0; i < N; i++ {
for j = 0; j < N; j++ {
ATA[i][j] = ATA[i][j] + H[i]*H[j]
}
ATb[i] = ATb[i] + H[i]*oc
}
}
func LeastSquareEquation(ATA [][]float64, ATb []float64) (X []float64) {
/*------------------------------------------------------
Created : Song Yezhi 2022-5-21 20:04
solve Least Square Equation by LDL'
--------------------------------------------------------
Input Parameters :
Output Parameters :
--------------------------------------------------------
Email : song.yz@foxmail.com
Copyrigt (C) : Chinese Academy of Sciences
All rights reserved, 2022
-------------------------------------------------------*/
N := len(ATA)
L := Mat(N, N)
D := Vec(N)
L, D = LDL(ATA)
X = Vec(N)
var tmp1 float64
var i, k int
y := Vec(N)
y[0] = ATb[0]
for i = 1; i < N; i++ {
tmp1 = 0.0
for k = 0; k <= i-1; k++ {
tmp1 = tmp1 + L[i][k]*y[k]
}
y[i] = ATb[i] - tmp1
}
X[N-1] = y[N-1] / D[N-1]
//-------------------------
//i = N - 1
for i=N-1; i>=0 ;i--{
tmp1 = 0.0
for k = i + 1; k < N; k++ {
tmp1 = tmp1 + L[k][i]*X[k]
}
X[i] = y[i]/D[i] - tmp1
/*----------------------
if i == 0 {
break
}
i = i-1
*/
}
return
}