1 文本格式
using System;
namespace Legalsoft.Truffer
{
/// <summary>
/// Object for interpolating missing data in a matrix by solving Laplace's
/// equation.Call constructor once, then solve one or more times
/// </summary>
public class Laplace_interp : Linbcg
{
private double[,] mat { get; set; }
private int ii { get; set; }
private int jj { get; set; }
private int nn { get; set; }
private int iter;
private double[] b;
private double[] y;
private double[] mask { get; set; }
/// <summary>
/// Values greater than 1.e99 in the input matrix mat are deemed to be missing
/// data.The matrix is not altered until solve is called.
/// </summary>
/// <param name="matrix"></param>
public Laplace_interp(double[,] matrix)
{
this.mat = matrix;
this.ii = mat.GetLength(0);
this.jj = mat.GetLength(1);
this.nn = ii * jj;
this.iter = 0;
this.b = new double[nn];
this.y = new double[nn];
this.mask = new double[nn];
double vl = 0.0;
for (int k = 0; k < nn; k++)
{
int i = k / jj;
int j = k - i * jj;
if (mat[i, j] < 1.0e99)
{
b[k] = y[k] = vl = mat[i, j];
mask[k] = 1;
}
else
{
b[k] = 0.0;
y[k] = vl;
mask[k] = 0;
}
}
}
/// <summary>
/// Diagonal preconditioner. (Diagonal elements all unity.)
/// </summary>
/// <param name="b"></param>
/// <param name="x"></param>
/// <param name="itrnsp"></param>
public override void asolve(double[] b, double[] x, int itrnsp)
{
int n = b.Length;
for (int i = 0; i < n; i++)
{
x[i] = b[i];
}
}
/// <summary>
/// Sparse matrix, and matrix transpose, multiply.
/// </summary>
/// <param name="x"></param>
/// <param name="r"></param>
/// <param name="itrnsp"></param>
public override void atimes(double[] x, double[] r, int itrnsp)
{
int n = r.Length;
double del;
for (int k = 0; k < n; k++)
{
r[k] = 0.0;
}
for (int k = 0; k < n; k++)
{
int i = k / jj;
int j = k - i * jj;
if (mask[k] > 0.0)
{
r[k] += x[k];
}
else if (i > 0 && i < ii - 1 && j > 0 && j < jj - 1)
{
if (itrnsp != 0)
{
r[k] += x[k];
del = -0.25 * x[k];
r[k - 1] += del;
r[k + 1] += del;
r[k - jj] += del;
r[k + jj] += del;
}
else
{
r[k] = x[k] - 0.25 * (x[k - 1] + x[k + 1] + x[k + jj] + x[k - jj]);
}
}
else if (i > 0 && i < ii - 1)
{
if (itrnsp != 0)
{
r[k] += x[k];
del = -0.5 * x[k];
r[k - jj] += del;
r[k + jj] += del;
}
else
{
r[k] = x[k] - 0.5 * (x[k + jj] + x[k - jj]);
}
}
else if (j > 0 && j < jj - 1)
{
if (itrnsp != 0)
{
r[k] += x[k];
del = -0.5 * x[k];
r[k - 1] += del;
r[k + 1] += del;
}
else
{
r[k] = x[k] - 0.5 * (x[k + 1] + x[k - 1]);
}
}
else
{
int jjt = i == 0 ? jj : -jj;
int it = j == 0 ? 1 : -1;
if (itrnsp != 0)
{
r[k] += x[k];
del = -0.5 * x[k];
r[k + jjt] += del;
r[k + it] += del;
}
else
{
r[k] = x[k] - 0.5 * (x[k + jjt] + x[k + it]);
}
}
}
}
/// <summary>
/// Invoke Linbcg::solve with appropriate arguments.The default argument
/// values will usually work, in which case this routine need be called only
/// once. The original matrix mat is refilled with the interpolated solution.
/// </summary>
/// <param name="tol"></param>
/// <param name="itmax"></param>
/// <returns></returns>
public double solve(double tol = 1.0e-6, int itmax = -1)
{
double err = 0.0;
if (itmax <= 0)
{
itmax = 2 * Math.Max(ii, jj);
}
solve( b, y, 1, tol, itmax, ref iter, ref err);
for (int k = 0, i = 0; i < ii; i++)
{
for (int j = 0; j < jj; j++)
{
mat[i, j] = y[k++];
}
}
return err;
}
}
}
2 代码格式
using System;namespace Legalsoft.Truffer
{/// <summary>/// Object for interpolating missing data in a matrix by solving Laplace's/// equation.Call constructor once, then solve one or more times/// </summary>public class Laplace_interp : Linbcg{private double[,] mat { get; set; }private int ii { get; set; }private int jj { get; set; }private int nn { get; set; }private int iter;private double[] b;private double[] y;private double[] mask { get; set; }/// <summary>/// Values greater than 1.e99 in the input matrix mat are deemed to be missing/// data.The matrix is not altered until solve is called./// </summary>/// <param name="matrix"></param>public Laplace_interp(double[,] matrix){this.mat = matrix;this.ii = mat.GetLength(0);this.jj = mat.GetLength(1);this.nn = ii * jj;this.iter = 0;this.b = new double[nn];this.y = new double[nn];this.mask = new double[nn];double vl = 0.0;for (int k = 0; k < nn; k++){int i = k / jj;int j = k - i * jj;if (mat[i, j] < 1.0e99){b[k] = y[k] = vl = mat[i, j];mask[k] = 1;}else{b[k] = 0.0;y[k] = vl;mask[k] = 0;}}}/// <summary>/// Diagonal preconditioner. (Diagonal elements all unity.)/// </summary>/// <param name="b"></param>/// <param name="x"></param>/// <param name="itrnsp"></param>public override void asolve(double[] b, double[] x, int itrnsp){int n = b.Length;for (int i = 0; i < n; i++){x[i] = b[i];}}/// <summary>/// Sparse matrix, and matrix transpose, multiply./// </summary>/// <param name="x"></param>/// <param name="r"></param>/// <param name="itrnsp"></param>public override void atimes(double[] x, double[] r, int itrnsp){int n = r.Length;double del;for (int k = 0; k < n; k++){r[k] = 0.0;}for (int k = 0; k < n; k++){int i = k / jj;int j = k - i * jj;if (mask[k] > 0.0){r[k] += x[k];}else if (i > 0 && i < ii - 1 && j > 0 && j < jj - 1){if (itrnsp != 0){r[k] += x[k];del = -0.25 * x[k];r[k - 1] += del;r[k + 1] += del;r[k - jj] += del;r[k + jj] += del;}else{r[k] = x[k] - 0.25 * (x[k - 1] + x[k + 1] + x[k + jj] + x[k - jj]);}}else if (i > 0 && i < ii - 1){if (itrnsp != 0){r[k] += x[k];del = -0.5 * x[k];r[k - jj] += del;r[k + jj] += del;}else{r[k] = x[k] - 0.5 * (x[k + jj] + x[k - jj]);}}else if (j > 0 && j < jj - 1){if (itrnsp != 0){r[k] += x[k];del = -0.5 * x[k];r[k - 1] += del;r[k + 1] += del;}else{r[k] = x[k] - 0.5 * (x[k + 1] + x[k - 1]);}}else{int jjt = i == 0 ? jj : -jj;int it = j == 0 ? 1 : -1;if (itrnsp != 0){r[k] += x[k];del = -0.5 * x[k];r[k + jjt] += del;r[k + it] += del;}else{r[k] = x[k] - 0.5 * (x[k + jjt] + x[k + it]);}}}}/// <summary>/// Invoke Linbcg::solve with appropriate arguments.The default argument/// values will usually work, in which case this routine need be called only/// once. The original matrix mat is refilled with the interpolated solution./// </summary>/// <param name="tol"></param>/// <param name="itmax"></param>/// <returns></returns>public double solve(double tol = 1.0e-6, int itmax = -1){double err = 0.0;if (itmax <= 0){itmax = 2 * Math.Max(ii, jj);}solve( b, y, 1, tol, itmax, ref iter, ref err);for (int k = 0, i = 0; i < ii; i++){for (int j = 0; j < jj; j++){mat[i, j] = y[k++];}}return err;}}
}
