song.yz@foxmail.com wechat: math-box

计算数学

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
}