Golang Triplet.Put方法代码示例

// AddToKb adds element K to global Jacobian matrix Kb
func (o *ElastRod) AddToKb(Kb *la.Triplet, sol *Solution, firstIt bool) (err error) {
	for i, I := range o.Umap {
		for j, J := range o.Umap {
			Kb.Put(I, J, o.K[i][j])
		}
	}
	return
}

/*  Jacobian
    ========
        Calculates (with N=n-1):
            df0dx0, df0dx1, df0dx2, ... df0dxN
            df1dx0, df1dx1, df1dx2, ... df1dxN
                 . . . . . . . . . . . . .
            dfNdx0, dfNdx1, dfNdx2, ... dfNdxN
    INPUT:
        ffcn : f(x) function
        x    : station where dfdx has to be calculated
        fx   : f @ x
        w    : workspace with size == n == len(x)
    RETURNS:
        J : dfdx @ x [must be pre-allocated]        */
func Jacobian(J *la.Triplet, ffcn Cb_f, x, fx, w []float64, distr bool) (err error) {
	ndim := len(x)
	start, endp1 := 0, ndim
	if distr {
		id, sz := mpi.Rank(), mpi.Size()
		start, endp1 = (id*ndim)/sz, ((id+1)*ndim)/sz
		if J.Max() == 0 {
			J.Init(ndim, ndim, (endp1-start)*ndim)
		}
	} else {
		if J.Max() == 0 {
			J.Init(ndim, ndim, ndim*ndim)
		}
	}
	J.Start()
	// NOTE: cannot split calculation by columns unless the f function is
	//       independently calculated by each MPI processor.
	//       Otherwise, the AllReduce in f calculation would
	//       join pieces of f from different processors calculated for
	//       different x values (δx[col] from different columns).
	/*
	   for col := start; col < endp1; col++ {
	       xsafe := x[col]
	       delta := math.Sqrt(EPS * max(CTE1, math.Abs(xsafe)))
	       x[col] = xsafe + delta
	       ffcn(w, x) // fnew
	       io.Pforan("x = %v, f = %v\n", x, w)
	       for row := 0; row < ndim; row++ {
	           J.Put(row, col, (w[row]-fx[row])/delta)
	       }
	       x[col] = xsafe
	   }
	*/
	var df float64
	for col := 0; col < ndim; col++ {
		xsafe := x[col]
		delta := math.Sqrt(EPS * max(CTE1, math.Abs(xsafe)))
		x[col] = xsafe + delta
		err = ffcn(w, x) // w := f(x+δx[col])
		if err != nil {
			return
		}
		for row := start; row < endp1; row++ {
			df = w[row] - fx[row]
			//if math.Abs(df) > EPS {
			J.Put(row, col, df/delta)
			//}
		}
		x[col] = xsafe
	}
	return
}

// AddToKb adds element K to global Jacobian matrix Kb
func (o *Rod) AddToKb(Kb *la.Triplet, sol *Solution, firstIt bool) (err error) {

	// zero K matrix
	la.MatFill(o.K, 0)
	la.MatFill(o.M, 0) // TODO: implement mass matrix

	// for each integration point
	var E float64
	nverts := o.Cell.Shp.Nverts
	for idx, ip := range o.IpsElem {

		// interpolation functions, gradients and variables @ ip
		err = o.ipvars(idx, sol)
		if err != nil {
			return
		}

		// auxiliary
		coef := ip[3]
		Jvec := o.Cell.Shp.Jvec3d
		G := o.Cell.Shp.Gvec
		J := o.Cell.Shp.J

		// add contribution to consistent tangent matrix
		for m := 0; m < nverts; m++ {
			for n := 0; n < nverts; n++ {
				for i := 0; i < o.Ndim; i++ {
					for j := 0; j < o.Ndim; j++ {
						r := i + m*o.Ndim
						c := j + n*o.Ndim
						E, err = o.Model.CalcD(o.States[idx], firstIt)
						if err != nil {
							return
						}
						o.K[r][c] += coef * o.A * E * G[m] * G[n] * Jvec[i] * Jvec[j] / J
					}
				}
			}
		}
	}

	// add K to sparse matrix Kb
	for i, I := range o.Umap {
		for j, J := range o.Umap {
			Kb.Put(I, J, o.K[i][j])
		}
	}
	return
}

func main() {

	mpi.Start(false)
	defer func() {
		mpi.Stop(false)
	}()

	if mpi.Rank() == 0 {
		chk.PrintTitle("Test SumToRoot 01")
	}

	M := [][]float64{
		{1000, 1000, 1000, 1011, 1021, 1000},
		{1000, 1000, 1000, 1012, 1022, 1000},
		{1000, 1000, 1000, 1013, 1023, 1000},
		{1011, 1012, 1013, 1000, 1000, 1000},
		{1021, 1022, 1023, 1000, 1000, 1000},
		{1000, 1000, 1000, 1000, 1000, 1000},
	}

	id, sz, m := mpi.Rank(), mpi.Size(), len(M)
	start, endp1 := (id*m)/sz, ((id+1)*m)/sz

	if sz > 6 {
		chk.Panic("this test works with at most 6 processors")
	}

	var J la.Triplet
	J.Init(m, m, m*m)
	for i := start; i < endp1; i++ {
		for j := 0; j < m; j++ {
			J.Put(i, j, M[i][j])
		}
	}
	la.PrintMat(fmt.Sprintf("J @ proc # %d", id), J.ToMatrix(nil).ToDense(), "%10.1f", false)

	la.SpTriSumToRoot(&J)
	var tst testing.T
	if mpi.Rank() == 0 {
		chk.Matrix(&tst, "J @ proc 0", 1.0e-17, J.ToMatrix(nil).ToDense(), [][]float64{
			{1000, 1000, 1000, 1011, 1021, 1000},
			{1000, 1000, 1000, 1012, 1022, 1000},
			{1000, 1000, 1000, 1013, 1023, 1000},
			{1011, 1012, 1013, 1000, 1000, 1000},
			{1021, 1022, 1023, 1000, 1000, 1000},
			{1000, 1000, 1000, 1000, 1000, 1000},
		})
	}
}

// adds element K to global Jacobian matrix Kb
func (o Rod) AddToKb(Kb *la.Triplet, sol *Solution, firstIt bool) (ok bool) {

	// zero K matrix
	la.MatFill(o.K, 0)
	la.MatFill(o.M, 0)

	// for each integration point
	nverts := o.Shp.Nverts
	ndim := Global.Ndim
	for idx, ip := range o.IpsElem {

		// interpolation functions, gradients and variables @ ip
		if !o.ipvars(idx, sol) {
			return
		}

		// auxiliary
		coef := ip.W
		Jvec := o.Shp.Jvec3d
		G := o.Shp.Gvec
		J := o.Shp.J

		// add contribution to consistent tangent matrix
		for m := 0; m < nverts; m++ {
			for n := 0; n < nverts; n++ {
				for i := 0; i < ndim; i++ {
					for j := 0; j < ndim; j++ {
						r := i + m*ndim
						c := j + n*ndim
						E, err := o.Model.CalcD(o.States[idx], firstIt)
						if LogErr(err, "AddToKb") {
							return
						}
						o.K[r][c] += coef * o.A * E * G[m] * G[n] * Jvec[i] * Jvec[j] / J
					}
				}
			}
		}
	}

	// add K to sparse matrix Kb
	for i, I := range o.Umap {
		for j, J := range o.Umap {
			Kb.Put(I, J, o.K[i][j])
		}
	}
	return true
}

// adds element K to global Jacobian matrix Kb
func (o Beam) AddToKb(Kb *la.Triplet, sol *Solution, firstIt bool) (ok bool) {
	if Global.Sim.Data.Steady {
		for i, I := range o.Umap {
			for j, J := range o.Umap {
				Kb.Put(I, J, o.K[i][j])
			}
		}
	} else {
		dc := Global.DynCoefs
		for i, I := range o.Umap {
			for j, J := range o.Umap {
				Kb.Put(I, J, o.M[i][j]*dc.α1+o.K[i][j])
			}
		}
	}
	return true
}

// adds element K to global Jacobian matrix Kb
func (o *Beam) AddToKb(Kb *la.Triplet, sol *Solution, firstIt bool) (err error) {
	if sol.Steady {
		for i, I := range o.Umap {
			for j, J := range o.Umap {
				Kb.Put(I, J, o.K[i][j])
			}
		}
		return
	}
	α1 := sol.DynCfs.α1
	for i, I := range o.Umap {
		for j, J := range o.Umap {
			Kb.Put(I, J, o.M[i][j]*α1+o.K[i][j])
		}
	}
	return
}

func main() {

	mpi.Start(false)
	defer func() {
		mpi.Stop(false)
	}()

	myrank := mpi.Rank()
	if myrank == 0 {
		chk.PrintTitle("Test MUMPS Sol 02")
	}

	ndim := 10
	id, sz := mpi.Rank(), mpi.Size()
	start, endp1 := (id*ndim)/sz, ((id+1)*ndim)/sz

	if mpi.Size() > ndim {
		chk.Panic("the number of processors must be smaller than or equal to %d", ndim)
	}

	b := make([]float64, ndim)
	var t la.Triplet
	t.Init(ndim, ndim, ndim*ndim)

	for i := start; i < endp1; i++ {
		j := i
		if i > 0 {
			j = i - 1
		}
		for ; j < ndim; j++ {
			val := 10.0 - float64(j)
			if i > j {
				val -= 1.0
			}
			t.Put(i, j, val)
		}
		b[i] = float64(i + 1)
	}

	x_correct := []float64{-1, 8, -65, 454, -2725, 13624, -54497, 163490, -326981, 326991}
	sum_b_to_root := true
	la.RunMumpsTestR(&t, 1e-4, b, x_correct, sum_b_to_root)
}

// AddToKb adds element K to global Jacobian matrix Kb
func (o *ElemPhi) AddToKb(Kb *la.Triplet, sol *Solution, firstIt bool) (err error) {

	// auxiliary
	nverts := o.Cell.Shp.Nverts

	// zero K matrix
	la.MatFill(o.K, 0)
	dt := sol.Dt * 2

	// for each integration point
	for iip, ip := range o.IpsElem {

		// interpolation functions and gradients
		err = o.Cell.Shp.CalcAtIp(o.X, ip, true)
		if err != nil {
			return
		}

		// auxiliary variables
		coef := o.Cell.Shp.J * ip[3]
		S := o.Cell.Shp.S
		G := o.Cell.Shp.G

		// add to right hand side vector
		for m := 0; m < nverts; m++ {
			for n := 0; n < nverts; n++ {
				o.K[m][n] -= coef * (S[m]*S[n] + dt*dt/6.0*(o.v_0[iip][0]*G[m][0]+o.v_0[iip][1]*G[m][1])*(o.v_0[iip][0]*G[n][0]+o.v_0[iip][1]*G[n][1])) / dt
			}
		}
	}

	// add K to sparse matrix Kb
	for i, I := range o.Umap {
		for j, J := range o.Umap {
			Kb.Put(I, J, o.K[i][j])
		}
	}
	return
}

func Assemble(K11, K12 *la.Triplet, F1 []float64, src Cb_src, g *Grid2D, e *Equations) {
	K11.Start()
	K12.Start()
	la.VecFill(F1, 0.0)
	kx, ky := 1.0, 1.0
	alp, bet, gam := 2.0*(kx/g.Dxx+ky/g.Dyy), -kx/g.Dxx, -ky/g.Dyy
	mol := []float64{alp, bet, bet, gam, gam}
	for i, I := range e.RF1 {
		col, row := I%g.Nx, I/g.Nx
		nodes := []int{I, I - 1, I + 1, I - g.Nx, I + g.Nx} // I, left, right, bottom, top
		if col == 0 {
			nodes[1] = nodes[2]
		}
		if col == g.Nx-1 {
			nodes[2] = nodes[1]
		}
		if row == 0 {
			nodes[3] = nodes[4]
		}
		if row == g.Ny-1 {
			nodes[4] = nodes[3]
		}
		for k, J := range nodes {
			j1, j2 := e.FR1[J], e.FR2[J] // 1 or 2?
			if j1 > -1 {                 // 11
				K11.Put(i, j1, mol[k])
			} else { // 12
				K12.Put(i, j2, mol[k])
			}
		}
		if src != nil {
			x := float64(col) * g.Dx
			y := float64(row) * g.Dy
			F1[i] += src(x, y)
		}
	}
}

// AddToKb adds element K to global Jacobian matrix Kb
func (o *ElemPhi) AddToKb(Kb *la.Triplet, sol *Solution, firstIt bool) (ok bool) {

	// auxiliary
	β1 := Global.DynCoefs.β1
	nverts := o.Shp.Nverts

	// zero K matrix
	la.MatFill(o.K, 0)

	// for each integration point
	for _, ip := range o.IpsElem {

		// interpolation functions and gradients
		if LogErr(o.Shp.CalcAtIp(o.X, ip, true), "InterpStarVars") {
			return
		}

		// auxiliary variables
		coef := o.Shp.J * ip.W
		S := o.Shp.S

		// add to right hand side vector
		for m := 0; m < nverts; m++ {
			for n := 0; n < nverts; n++ {
				o.K[m][n] += coef * S[m] * S[n] * β1
			}
		}
	}

	// add K to sparse matrix Kb
	for i, I := range o.Umap {
		for j, J := range o.Umap {
			Kb.Put(I, J, o.K[i][j])
		}
	}
	return true
}

// Jacobian computes Jacobian (sparse) matrix
//      Calculates (with N=n-1):
//          df0dx0, df0dx1, df0dx2, ... df0dxN
//          df1dx0, df1dx1, df1dx2, ... df1dxN
//               . . . . . . . . . . . . .
//          dfNdx0, dfNdx1, dfNdx2, ... dfNdxN
//  INPUT:
//      ffcn : f(x) function
//      x    : station where dfdx has to be calculated
//      fx   : f @ x
//      w    : workspace with size == n == len(x)
//  RETURNS:
//      J : dfdx @ x [must be pre-allocated]
func Jacobian(J *la.Triplet, ffcn Cb_f, x, fx, w []float64) (err error) {
	ndim := len(x)
	start, endp1 := 0, ndim
	if J.Max() == 0 {
		J.Init(ndim, ndim, ndim*ndim)
	}
	J.Start()
	var df float64
	for col := 0; col < ndim; col++ {
		xsafe := x[col]
		delta := math.Sqrt(EPS * max(CTE1, math.Abs(xsafe)))
		x[col] = xsafe + delta
		err = ffcn(w, x) // w := f(x+δx[col])
		if err != nil {
			return
		}
		for row := start; row < endp1; row++ {
			df = w[row] - fx[row]
			J.Put(row, col, df/delta)
		}
		x[col] = xsafe
	}
	return
}

func IpBmatrix_sparse(B *la.Triplet, ndim, nne int, G [][]float64, radius float64, S []float64, axisym bool) {
	B.Start()
	if ndim == 3 {
		for i := 0; i < nne; i++ {
			B.Put(0, 0+i*3, G[i][0])
			B.Put(1, 1+i*3, G[i][1])
			B.Put(2, 2+i*3, G[i][2])
			B.Put(3, 0+i*3, G[i][1]/SQ2)
			B.Put(4, 1+i*3, G[i][2]/SQ2)
			B.Put(5, 2+i*3, G[i][0]/SQ2)
			B.Put(3, 1+i*3, G[i][0]/SQ2)
			B.Put(4, 2+i*3, G[i][1]/SQ2)
			B.Put(5, 0+i*3, G[i][2]/SQ2)
		}
		return
	}
	if axisym {
		for i := 0; i < nne; i++ {
			B.Put(0, 0+i*2, G[i][0])
			B.Put(1, 1+i*2, G[i][1])
			B.Put(2, 0+i*2, S[i]/radius)
			B.Put(3, 0+i*2, G[i][1]/SQ2)
			B.Put(3, 1+i*2, G[i][0]/SQ2)
		}
		return
	}
	for i := 0; i < nne; i++ {
		B.Put(0, 0+i*2, G[i][0])
		B.Put(1, 1+i*2, G[i][1])
		B.Put(3, 0+i*2, G[i][1]/SQ2)
		B.Put(3, 1+i*2, G[i][0]/SQ2)
	}
}

func main() {

	// input matrix in Triplet format
	// including repeated positions. e.g. (0,0)
	var A la.Triplet
	A.Init(5, 5, 13)
	A.Put(0, 0, 1.0) // << repeated
	A.Put(0, 0, 1.0) // << repeated
	A.Put(1, 0, 3.0)
	A.Put(0, 1, 3.0)
	A.Put(2, 1, -1.0)
	A.Put(4, 1, 4.0)
	A.Put(1, 2, 4.0)
	A.Put(2, 2, -3.0)
	A.Put(3, 2, 1.0)
	A.Put(4, 2, 2.0)
	A.Put(2, 3, 2.0)
	A.Put(1, 4, 6.0)
	A.Put(4, 4, 1.0)

	// right-hand-side
	b := []float64{8.0, 45.0, -3.0, 3.0, 19.0}

	// solve
	x, err := la.SolveRealLinSys(&A, b)
	if err != nil {
		io.Pfred("solver failed:\n%v", err)
		return
	}

	// output
	la.PrintMat("a", A.ToMatrix(nil).ToDense(), "%5g", false)
	la.PrintVec("b", b, "%v ", false)
	la.PrintVec("x", x, "%v ", false)
}

func main() {

	// input matrix in Triplet format
	// including repeated positions. e.g. (0,0)
	var A la.Triplet
	A.Init(5, 5, 13)
	A.Put(0, 0, 1.0) // << repeated
	A.Put(0, 0, 1.0) // << repeated
	A.Put(1, 0, 3.0)
	A.Put(0, 1, 3.0)
	A.Put(2, 1, -1.0)
	A.Put(4, 1, 4.0)
	A.Put(1, 2, 4.0)
	A.Put(2, 2, -3.0)
	A.Put(3, 2, 1.0)
	A.Put(4, 2, 2.0)
	A.Put(2, 3, 2.0)
	A.Put(1, 4, 6.0)
	A.Put(4, 4, 1.0)

	// right-hand-side
	b := []float64{8.0, 45.0, -3.0, 3.0, 19.0}

	// allocate solver
	lis := la.GetSolver("umfpack")
	defer lis.Clean()

	// info
	symmetric := false
	verbose := false
	timing := false

	// initialise solver (R)eal
	err := lis.InitR(&A, symmetric, verbose, timing)
	if err != nil {
		io.Pfred("solver failed:\n%v", err)
		return
	}

	// factorise
	err = lis.Fact()
	if err != nil {
		io.Pfred("solver failed:\n%v", err)
		return
	}

	// solve (R)eal
	var dummy bool
	x := make([]float64, len(b))
	err = lis.SolveR(x, b, dummy) // x := inv(a) * b
	if err != nil {
		io.Pfred("solver failed:\n%v", err)
		return
	}

	// output
	la.PrintMat("a", A.ToMatrix(nil).ToDense(), "%5g", false)
	la.PrintVec("b", b, "%v ", false)
	la.PrintVec("x", x, "%v ", false)
}