using System; using System.Collections.Generic; using System.ComponentModel; using System.Data; using System.Drawing; using System.Linq; using System.Text; using System.Windows.Forms; namespace vcsGJACOBI { public partial class Form1 : Form { public Form1() { InitializeComponent(); } //------------------------------------------------------------------ private void toolStripButton1_Click(object sender, EventArgs e) { int i, j, k, nf; double[,] a; double[,] b; double[] ev; double[,] vec; string dat; double[,] ar; double[,] br; double[,] c; double[,] d; int jdg; nf = 3; a = new double[nf + 1, nf + 1]; b = new double[nf + 1, nf + 1]; ev = new double[nf + 1]; vec = new double[nf + 1, nf + 1]; ar = new double[nf + 1, nf + 1]; br = new double[nf + 1, nf + 1]; c = new double[nf + 1, nf + 1]; d = new double[nf + 1, nf + 1]; //Case-1 //a[1, 1] = 4; a[1, 2] = 2; a[1, 3] = 3; //a[2, 1] = 2; a[2, 2] = 6; a[2, 3] = 1; //a[3, 1] = 3; a[3, 2] = 1; a[3, 3] = 5; //b[1, 1] = 1; b[1, 2] = 0; b[1, 3] = 0; //b[2, 1] = 0; b[2, 2] = 1; b[2, 3] = 0; //b[3, 1] = 0; b[3, 2] = 0; b[3, 3] = 1; a[1, 1] = 106.482 ; a[1, 2] = -59.228 ; a[1, 3] = 10.816; a[2, 1] = -59.228 ; a[2, 2] = 89.1 ; a[2, 3] = -42.618; a[3, 1] = 10.816 ; a[3, 2] = -42.618 ; a[3, 3] = 33.348; b[1, 1] = 0.01 ; b[1, 2] = 0.0 ; b[1, 3] = 0.0 ; b[2, 1] = 0.0 ; b[2, 2] = 0.01 ; b[2, 3] = 0.0 ; b[3, 1] = 0.0 ; b[3, 2] = 0.0 ; b[3, 3] = 0.01; //a[1, 1] = 3 ; a[1, 2] = 2 ; a[1, 3] = -2; //a[2, 1] = 2 ; a[2, 2] = 3 ; a[2, 3] = -2; //a[3, 1] = -2 ; a[3, 2] = -2 ; a[3, 3] = 8; //b[1, 1] = 1 ; b[1, 2] = 0 ; b[1, 3] = 0; //b[2, 1] = 0 ; b[2, 2] = 1 ; b[2, 3] = 0; //b[3, 1] = 0 ; b[3, 2] = 0 ; b[3, 3] = 1; //a[1, 1] = 5 ; a[1, 2] = 4 ; a[1, 3] = 3 ; a[1, 4] = 2 ; a[1, 5] = 1; //a[2, 1] = 4 ; a[2, 2] = 4 ; a[2, 3] = 3 ; a[2, 4] = 2 ; a[2, 5] = 1; //a[3, 1] = 3 ; a[3, 2] = 3 ; a[3, 3] = 3 ; a[3, 4] = 2 ; a[3, 5] = 1; //a[4, 1] = 2 ; a[4, 2] = 2 ; a[4, 3] = 2 ; a[4, 4] = 2 ; a[4, 5] = 1; //a[5, 1] = 1 ; a[5, 2] = 1 ; a[5, 3] = 1 ; a[5, 4] = 1 ; a[5, 5] = 1; //b[1, 1] = 1 ; b[1, 2] = 0 ; b[1, 3] = 0 ; b[1, 4] = 0 ; b[1, 5] = 0; //b[2, 1] = 0 ; b[2, 2] = 1 ; b[2, 3] = 0 ; b[2, 4] = 0 ; b[2, 5] = 0; //b[3, 1] = 0 ; b[3, 2] = 0 ; b[3, 3] = 1 ; b[3, 4] = 0 ; b[3, 5] = 0; //b[4, 1] = 0 ; b[4, 2] = 0 ; b[4, 3] = 0 ; b[4, 4] = 1 ; b[4, 5] = 0; //b[5, 1] = 0 ; b[5, 2] = 0 ; b[5, 3] = 0 ; b[5, 4] = 0 ; b[5, 5] = 1; for (i = 1; i <= nf; i++) { for (j = 1; j <= nf; j++) { ar[i, j] = a[i, j]; br[i, j] = b[i, j]; } } GJACOBI(nf, a, b, ev, vec); Console.WriteLine("*** Eigen value ***"); dat = String.Format("{0,16:0.000000e+000}", ev[1]); for (j = 2; j <= nf; j++) { dat = dat + String.Format("{0,16:0.000000e+000}", ev[j]); } Console.WriteLine(dat); Console.WriteLine("*** Eigen vector ***"); for (i = 1; i <= nf; i++) { dat = String.Format("{0,16:0.000000e+000}", vec[i, 1]); for (j = 2; j <= nf; j++) { dat = dat + String.Format("{0,16:0.000000e+000}", vec[i, j]); } Console.WriteLine(dat); } //Case-1 出力結果 //*** Eigen value *** //9.000000e+000 4.732051e+000 1.267949e+000 //*** Eigen vector *** //9.227825e-001 -2.633891e-001 1.060632e+000 //9.227825e-001 9.829816e-001 -2.841954e-001 //9.227825e-001 -7.195925e-001 -7.764364e-001 //検算 for (i = 1; i <= nf; i++) { for (j = 1; j <= nf; j++) { c[i, j] = 0.0; d[i, j] = 0.0; for (k = 1; k <= nf; k++) { c[i, j] = c[i, j] + ar[i, k] * vec[k, j]; d[i, j] = d[i, j] + br[i, k] * vec[k, j] * ev[j]; } } } Console.WriteLine("*** 検算 ***"); for (i = 1; i <= nf; i++) { dat = String.Format("{0,16:0.000000e+000}", c[i, 1]); for (j = 2; j <= nf; j++) { dat = dat + String.Format("{0,16:0.000000e+000}", c[i, j]); } dat = dat + " "; dat = dat + String.Format("{0,16:0.000000e+000}", d[i, 1]); for (j = 2; j <= nf; j++) { dat = dat + String.Format("{0,16:0.000000e+000}", d[i, j]); } Console.WriteLine(dat); } jdg = GJACOBIF(nf, a, b, ev, vec); Console.WriteLine("*** Eigen value ***"); dat = String.Format("{0,16:0.000000e+000}", ev[1]); for (j = 2; j <= nf; j++) { dat = dat + String.Format("{0,16:0.000000e+000}", ev[j]); } Console.WriteLine(dat); Console.WriteLine("*** Eigen vector ***"); for (i = 1; i <= nf; i++) { dat = String.Format("{0,16:0.000000e+000}", vec[i, 1]); for (j = 2; j <= nf; j++) { dat = dat + String.Format("{0,16:0.000000e+000}", vec[i, j]); } Console.WriteLine(dat); } //検算 for (i = 1; i <= nf; i++) { for (j = 1; j <= nf; j++) { c[i, j] = 0.0; d[i, j] = 0.0; for (k = 1; k <= nf; k++) { c[i, j] = c[i, j] + ar[i, k] * vec[k,j]; d[i, j] = d[i, j] + br[i, k] * vec[k,j] * ev[j]; } } } Console.WriteLine("*** 検算 ***"); for (i = 1; i <= nf; i++) { dat = String.Format("{0,16:0.000000e+000}", c[i, 1]); for (j = 2; j <= nf; j++) { dat = dat + String.Format("{0,16:0.000000e+000}", c[i, j]); } dat = dat + " "; dat = dat + String.Format("{0,16:0.000000e+000}", d[i, 1]); for (j = 2; j <= nf; j++) { dat = dat + String.Format("{0,16:0.000000e+000}", d[i, j]); } Console.WriteLine(dat); } } //------------------------------------------------------------------ private void GJACOBI(int nf, double[,] a, double[,] b, double[] ev, double[,] vec) { //************************************ //一般化JACOBI法による固有値計算 //[パソコン構造力学，宮澤健二著] //************************************ int i, j, k, kmax, ip, iq, im; double fmk, fmg, epsk, epsg; double fact = 1.0E-20; double amax, dam, de1, de2, de3, dd, ss, x, gg, ww; fmk = 0.0; fmg = 0.0; for (i = 1; i <= nf; i++) { for (j = 1; j <= nf; j++) { if (Math.Abs(a[i, j]) > fmk) fmk = Math.Abs(a[i, j]); if (Math.Abs(b[i, j]) > fmg) fmg = Math.Abs(b[i, j]); } } epsk = fmk * fact; epsg = fmg * fact; for (i = 1; i <= nf; i++) { for (j = 1; j <= nf; j++) { vec[i, j] = 0.0; } vec[i, i] = 1.0; } kmax = 5 * nf * nf; for (k = 1; k <= kmax; k++) { ip = 1; iq = 2; amax = a[1, 2] * a[1, 2] / (a[1, 1] * a[1, 1] + a[2, 2] * a[2, 2]); for (i = 1; i <= nf - 1; i++) { for (j = i + 1; j <= nf; j++) { dam = a[i, j] * a[i, j] / (a[i, i] * a[i, i] + a[j, j] * a[j, j]); if (amax < dam) { ip = i; iq = j; amax = dam; } } } for (i = 1; i <= nf - 1; i++) { for (j = i + 1; j <= nf; j++) { dam = b[i, j] * b[i, j] / (b[i, i] * b[i, i] + b[j, j] * b[j, j]); if (amax < dam) { ip = i; iq = j; amax = dam; } } } de1 = a[ip, ip] * b[ip, iq] - a[ip, iq] * b[ip, ip]; de2 = a[iq, iq] * b[ip, iq] - a[ip, iq] * b[iq, iq]; de3 = a[ip, ip] * b[iq, iq] - a[iq, iq] * b[ip, ip]; dd = de3 * de3 + 4.0 * de1 * de2; ss = Math.Sqrt(dd); if (de3 < 0) { x = 0.5 * (de3 - ss); } else { x = 0.5 * (de3 + ss); } if (x == 0.0) { dam = 0.0; gg = -a[ip, iq] / a[iq, iq]; } else { dam = de2 / x; gg = -de1 / x; } for (j = 1; j <= nf; j++) { ww = a[ip, j]; a[ip, j] = a[ip, j] + gg * a[iq, j]; a[iq, j] = a[iq, j] + dam * ww; ww = b[ip, j]; b[ip, j] = b[ip, j] + gg * b[iq, j]; b[iq, j] = b[iq, j] + dam * ww; } for (i = 1; i <= nf; i++) { ww = a[i, ip]; a[i, ip] = a[i, ip] + gg * a[i, iq]; a[i, iq] = a[i, iq] + dam * ww; ww = b[i, ip]; b[i, ip] = b[i, ip] + gg * b[i, iq]; b[i, iq] = b[i, iq] + dam * ww; ww = vec[i, ip]; vec[i, ip] = vec[i, ip] + gg * vec[i, iq]; vec[i, iq] = vec[i, iq] + dam * ww; } if ((amax < epsk) && (amax < epsg)) break; } for (i = 1; i <= nf; i++) { ev[i] = a[i, i] / b[i, i]; } //並び替え //ev[k]：最大固有値k=1，最小固有値k=nf //vec[i,k]：k番目の固有値に対応する固有ベクトル im = 0; for (k = 1; k <= nf - 1; k++) { amax = -1.0E+30; for (i = k; i <= nf; i++) { if (ev[i] >= amax) { im = i; amax = ev[i]; } } ww = ev[k]; ev[k] = ev[im]; ev[im] = ww; for (i = 1; i <= nf; i++) { ww = vec[i, k]; vec[i, k] = vec[i, im]; vec[i, im] = ww; } } } //------------------------------------------------------------------ private int GJACOBIF(int n, double[,] a, double[,] b, double[] eigv, double[,] x) { double[] d = new double[n + 1]; double rtol = 1e-30; int nsmax = 5*n*n; int i, j, k, jj; int nsweep, nr; double eps; double eptola, eptolb; double akk, ajj, ab, check; double sqch, d1, d2, den, ca, cg; int jp1, jm1, kp1, km1; double aj, bj, ak, bk; double xj, xk; double tol, dif; double epsa, epsb; double bb; //Initialize eigenvalue and eigenvector matrices for (i = 1; i <= n; i++) { if (a[i, i] < 0.0 || b[i, i] < 0.0) return 1; d[i] = a[i, i] / b[i, i]; eigv[i] = d[i]; } for (i = 1; i <= n; i++) { for (j = 1; j <= n; j++) { x[i, j] = 0.0; } x[i, i] = 1.0; } if (n == 1) return 0; //Initialize sweep counter and iteration nsweep = 0; nr = n - 1; num40: nsweep = nsweep + 1; //Check if present off-diagonal element is larger enough to require zeroing eps = Math.Pow(0.01, (double)nsweep) * Math.Pow(0.01, (double)nsweep); for (j = 1; j <= nr; j++) { jj = j + 1; for (k = jj; k <= n; k++) { eptola = (a[j, k] * a[j, k]) / (a[j, j] * a[k, k]); eptolb = (b[j, k] * b[j, k]) / (b[j, j] * b[k, k]); if (eptola < eps && eptolb < eps) continue; //If zeroing is required, calculate the rotation matrix elements ca and cg akk = a[k, k] * b[j, k] - b[k, k] * a[j, k]; ajj = a[j, j] * b[j, k] - b[j, j] * a[j, k]; ab = a[j, j] * b[k, k] - a[k, k] * b[j, j]; check = (ab * ab + 4.0 * akk * ajj) / 4.0; if (check < 0.0) return 1; sqch = Math.Sqrt(check); d1 = ab / 2.0 + sqch; d2 = ab / 2.0 - sqch; den = d1; if (Math.Abs(d1) < Math.Abs(d2)) den = d2; if (den == 0.0) { ca = 0.0; cg = -a[j, k] / a[k, k]; } else { ca = akk / den; cg = -ajj / den; } //Perform the generalized rotation to zero the present off-diagonal element if (n - 2 != 0) { jp1 = j + 1; jm1 = j - 1; kp1 = k + 1; km1 = k - 1; if (0 <= jm1 - 1) { for (i = 1; i <= jm1; i++) { aj = a[i, j]; bj = b[i, j]; ak = a[i, k]; bk = b[i, k]; a[i, j] = aj + cg * ak; b[i, j] = bj + cg * bk; a[i, k] = ak + ca * aj; b[i, k] = bk + ca * bj; } } if (kp1 - n <= 0) { for (i = kp1; i <= n; i++) { aj = a[j, i]; bj = b[j, i]; ak = a[k, i]; bk = b[k, i]; a[j, i] = aj + cg * ak; b[j, i] = bj + cg * bk; a[k, i] = ak + ca * aj; b[k, i] = bk + ca * bj; } } if (jp1 - km1 <= 0) { for (i = jp1; i <= km1; i++) { aj = a[j, i]; bj = b[j, i]; ak = a[i, k]; bk = b[i, k]; a[j, i] = aj + cg * ak; b[j, i] = bj + cg * bk; a[i, k] = ak + ca * aj; b[i, k] = bk + ca * bj; } } } ak = a[k, k]; bk = b[k, k]; a[k, k] = ak + 2.0 * ca * a[j, k] + ca * ca * a[j, j]; b[k, k] = bk + 2.0 * ca * b[j, k] + ca * ca * b[j, j]; a[j, j] = a[j, j] + 2.0 * cg * a[j, k] + cg * cg * ak; b[j, j] = b[j, j] + 2.0 * cg * b[j, k] + cg * cg * bk; a[j, k] = 0.0; b[j, k] = 0.0; //Update the eigenvector matrix after each rotation for (i = 1; i <= n; i++) { xj = x[i, j]; xk = x[i, k]; x[i, j] = xj + cg * xk; x[i, k] = xk + ca * xj; } } } //Update the eigenvalues after each sweep for (i = 1; i <= n; i++) { if (a[i, i] < 0.0 || b[i, i] < 0.0) return 1; eigv[i] = a[i, i] / b[i, i]; } //Check for convergence for (i = 1; i <= n; i++) { tol = rtol * d[i]; dif = Math.Abs(eigv[i] - d[i]); if (tol < dif) goto num280; } //Check all off-diagonal elements to see if another sweep is required eps = rtol * rtol; for (j = 1; j <= nr; j++) { jj = j + 1; for (k = jj; k <= n; k++) { epsa = (a[j, k] * a[j, k]) / (a[j, j] * a[k, k]); epsb = (b[j, k] * b[j, k]) / (b[j, j] * b[k, k]); if (eps < epsa || eps < epsb) goto num280; } } goto num255; //Update d vector and start new sweep, if allowed num280: for (i = 1; i <= n; i++) { d[i] = eigv[i]; } if (nsweep < nsmax) goto num40; //Fill out bottom triangle of resultant matrices and scale eigenvectors num255: for (i = 1; i <= n; i++) { for (j = 1; j <= n; j++) { a[j, i] = a[i, j]; b[j, i] = b[i, j]; } } for (j = 1; j <= n; j++) { bb = Math.Sqrt(b[j, j]); for (k = 1; k <= n; k++) { x[k, j] = x[k, j] / bb; } } return 0; } //------------------------------------------------------------------ } }