API tools faq
paste
Advertisement
Guest User

Untitled

a guest
Oct 16th, 2019
109
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 17.19 KB | None | 0 0
raw download clone embed print report
  1. package ru.omsu.imit.mgen;
  2.  
  3. import static java.lang.Math.abs;
  4.  
  5. /**
  6. * Hello world!
  7. */
  8. public class Gen {
  9. public Gen() {
  10. }
  11.  
  12.  
  13. public void print_matr(double[][] a, int n) {
  14. System.out.println(" ");
  15. System.out.println(" a:");
  16. for (int i = 0; i < n; i++) {
  17. for (int j = 0; j < n; j++) {
  18. System.out.print(a[i][j] + " ");
  19. }
  20. System.out.println();
  21. }
  22.  
  23. }
  24.  
  25.  
  26. public double matr_inf_norm(double[][] a, int n) {
  27.  
  28. int i, j;
  29. double s, norm = 0.;
  30.  
  31. for (i = 0; i < n; i++) {
  32. for (s = 0., j = 0; j < n; j++) s += Math.abs(a[i][j]);
  33. if (s > norm) norm = s;
  34. }
  35.  
  36. return norm;
  37. }
  38.  
  39.  
  40. public void matr_mul(double[][] a, double[][] b, double[][] c, int n) {
  41. int i, j, k;
  42.  
  43. for (i = 0; i < n; i++)
  44. for (j = 0; j < n; j++)
  45. for (c[i][j] = 0., k = 0; k < n; k++) c[i][j] += a[i][k] * b[k][j];
  46.  
  47. }
  48.  
  49.  
  50. public void Q_matrix(double[][] Q, int n, int schema) {
  51. int i, j;
  52. double q;
  53.  
  54. double curr, next = 1.;
  55. for (j = 0; j < n - 1; j++) {
  56. curr = next;
  57. next += 1.;
  58.  
  59. q = 1. / Math.sqrt(curr * next);
  60. for (i = 0; i <= j; i++) Q[i][j] = q;
  61. Q[j + 1][j] = -Math.sqrt(curr / next);
  62. for (i = j + 2; i < n; i++) Q[i][j] = 0.;
  63. }
  64.  
  65. q = 1. / Math.sqrt((double) n);
  66. for (i = 0; i < n; i++) Q[i][n - 1] = q;
  67. }
  68.  
  69.  
  70. public void mygen(double[][] a, double[][] a_inv, int n, double alpha, double beta, int sign_law, int lambda_law, int variant, int schema) {
  71. int i, j, k;
  72.  
  73. System.out.println(" M A T R I X G E N. ");
  74. System.out.println(" N = " + n);
  75. System.out.println(" | lambda_min | = " + alpha);
  76. System.out.println(" | lambda_max | = " + beta);
  77.  
  78. double[] lambda = new double[n];
  79.  
  80. // распределение знаков
  81. System.out.println(" sign_law = " + sign_law);
  82.  
  83. double[] sign = new double[n];
  84. for (i = 0; i < n; i++) sign[i] = 1.;
  85.  
  86. switch (sign_law) {
  87. case -1:
  88. for (i = 0; i < n; i++) sign[i] = -1.;
  89. break;
  90.  
  91. case 0:
  92. sign[0] = 1.;
  93. for (i = 1; i < n; i++) sign[i] = -sign[i - 1];
  94. break;
  95.  
  96. //другие законы распределения знаков
  97. // ...
  98.  
  99. }
  100. /* for( i=0; i<n; i++ ) cout<<sign[i]<<" ";
  101. cout<<endl;
  102. */
  103.  
  104. //распределение собственнных чисел
  105. System.out.println(" lambda_law = " + lambda_law);
  106.  
  107. double[] kappa = new double[n];
  108. for (i = 0; i < n; i++) kappa[i] = (double) i / (double) (n - 1);
  109. switch (lambda_law) {
  110. case 1:
  111. System.out.println(" kappa = sqrt( ) ");
  112. for (i = 0; i < n; i++) kappa[i] = Math.sqrt(kappa[i]);
  113. break;
  114.  
  115. case 2:
  116. System.out.println(" kappa = sin( ) ");
  117. double pi_half = Math.acos(-1.) * 0.5;
  118. for (i = 0; i < n; i++) kappa[i] = Math.sin(pi_half * kappa[i]);
  119. break;
  120.  
  121. //другие законы распределения собственных чисел
  122. // ...
  123.  
  124. }
  125. /* for( i=0; i<n; i++ ) cout<<kappa[i]<<" ";
  126. cout<<endl;
  127. */
  128.  
  129.  
  130. double[] J = new double[n];
  131. for (i = 0; i < n; i++) J[i] = sign[i] * ((1. - kappa[i]) * alpha + kappa[i] * beta);
  132. /* for( i=0; i<n; i++ ) cout<<J[i]<<" ";
  133. cout<<endl;
  134. */
  135.  
  136. double[] J_inv = new double[n];
  137. for (i = 0; i < n; i++) J_inv[i] = 1. / J[i];
  138.  
  139. double[][] Q = new double[n][];
  140. for (i = 0; i < n; i++) Q[i] = new double[n];
  141.  
  142. double[] aa = new double[3];
  143.  
  144.  
  145. System.out.println(" variant = " + variant);
  146.  
  147. switch (variant) {
  148. case 0: //симметричная матрица
  149. System.out.println(" simmetric matrix:");
  150. System.out.println(" schema = " + schema);
  151. switch (schema) {
  152. case 1:
  153. Q_matrix(Q, n, schema);
  154.  
  155. for (a[0][0] = 0., k = 0; k < n; k++) a[0][0] += Q[0][k] * J[k] * Q[0][k];
  156. for (j = 1; j < n; j++) {
  157. for (a[0][j] = 0., k = j - 1; k < n; k++) a[0][j] += Q[0][k] * J[k] * Q[j][k];
  158. a[j][0] = a[0][j];
  159. }
  160. for (i = 1; i < n; i++) {
  161. for (a[i][i] = 0., k = i - 1; k < n; k++) a[i][i] += Q[i][k] * J[k] * Q[i][k];
  162. for (j = i + 1; j < n; j++) {
  163. for (a[i][j] = 0., k = j - 1; k < n; k++) a[i][j] += Q[i][k] * J[k] * Q[j][k];
  164. a[j][i] = a[i][j];
  165. }
  166. }
  167.  
  168. //_______
  169. for (a_inv[0][0] = 0., k = 0; k < n; k++) a_inv[0][0] += Q[0][k] * J_inv[k] * Q[0][k];
  170. for (j = 1; j < n; j++) {
  171. for (a_inv[0][j] = 0., k = j - 1; k < n; k++) a_inv[0][j] += Q[0][k] * J_inv[k] * Q[j][k];
  172. a_inv[j][0] = a_inv[0][j];
  173. }
  174. for (i = 1; i < n; i++) {
  175. for (a_inv[i][i] = 0., k = i - 1; k < n; k++) a_inv[i][i] += Q[i][k] * J_inv[k] * Q[i][k];
  176. for (j = i + 1; j < n; j++) {
  177. for (a_inv[i][j] = 0., k = j - 1; k < n; k++)
  178. a_inv[i][j] += Q[i][k] * J_inv[k] * Q[j][k];
  179. a_inv[j][i] = a_inv[i][j];
  180. }
  181. }
  182. break;
  183.  
  184. }//schema
  185. break;
  186.  
  187. case 1: //матрица простой структуры
  188. System.out.println(" simple structure matrix:");
  189. System.out.println(" schema = " + schema);
  190. switch (schema) {
  191. case 1:
  192. //TJ
  193. //первая строка
  194. a[0][0] = J[0];
  195. a[0][1] = -J[1];
  196. // for(j=2; j<n; j++ ) a[0][j] = 0.;
  197. //до последней
  198. for (i = 1; i < n - 1; i++) {
  199. // for(j=0; j<i-1; j++ ) a[i][j] = 0.;
  200. a[i][i - 1] = -J[i - 1];
  201. a[i][i] = J[i] + J[i];
  202. a[i][i + 1] = -J[i + 1];
  203. // for(j=i+2; j<n; j++ ) a[i][j] = 0.;
  204. }
  205. //последняя (n-1)
  206. // for(j=0; j<n-2; j++ ) a[n-1][j] = 0.;
  207. a[n - 1][n - 2] = -J[n - 2];
  208. a[n - 1][n - 1] = J[n - 1] + J[n - 1];
  209.  
  210. //(TJ)T^{-1}
  211. //первая строка
  212. aa[1] = a[0][0];
  213. aa[2] = a[0][1];
  214. a[0][0] = aa[1] * (double) n + aa[2] * (double) (n - 1);
  215. double s = aa[1] + aa[2];
  216. for (j = 1; j < n; j++) a[0][j] = s * (double) (n - j);
  217. //до последней
  218. for (i = 1; i < n - 1; i++) {
  219. aa[0] = a[i][i - 1];
  220. aa[1] = a[i][i];
  221. aa[2] = a[i][i + 1];
  222. for (j = 0; j < i; j++)
  223. a[i][j] = aa[0] * (double) (n - i + 1) + aa[1] * (double) (n - i) + aa[2] * (double) (n - i - 1);
  224. s = aa[0] + aa[1];
  225. a[i][i] = s * (double) (n - i) + aa[2] * (double) (n - i - 1);
  226. s += aa[2];
  227. for (j = i + 1; j < n; j++) a[i][j] = s * (double) (n - j);
  228. }
  229. //последняя (n-1)
  230. aa[0] = a[n - 1][n - 2];
  231. aa[1] = a[n - 1][n - 1];
  232. s = aa[0] + aa[0] + aa[1];
  233. for (j = 0; j < n - 1; j++) a[n - 1][j] = s;
  234. a[n - 1][n - 1] = aa[0] + aa[1];
  235. //_______
  236.  
  237. //TJ^{-1}
  238. //первая строка
  239. a_inv[0][0] = J_inv[0];
  240. a_inv[0][1] = -J_inv[1];
  241. //до последней
  242. for (i = 1; i < n - 1; i++) {
  243. a_inv[i][i - 1] = -J_inv[i - 1];
  244. a_inv[i][i] = J_inv[i] + J_inv[i];
  245. a_inv[i][i + 1] = -J_inv[i + 1];
  246. }
  247. //последняя (n-1)
  248. a_inv[n - 1][n - 2] = -J_inv[n - 2];
  249. a_inv[n - 1][n - 1] = J_inv[n - 1] + J_inv[n - 1];
  250.  
  251. //(TJ^{-1})T^{-1}
  252. //первая строка
  253. aa[1] = a_inv[0][0];
  254. aa[2] = a_inv[0][1];
  255. a_inv[0][0] = aa[1] * (double) n + aa[2] * (double) (n - 1);
  256. s = aa[1] + aa[2];
  257. for (j = 1; j < n; j++) a_inv[0][j] = s * (double) (n - j);
  258. //до последней
  259. for (i = 1; i < n - 1; i++) {
  260. aa[0] = a_inv[i][i - 1];
  261. aa[1] = a_inv[i][i];
  262. aa[2] = a_inv[i][i + 1];
  263. for (j = 0; j < i; j++)
  264. a_inv[i][j] = aa[0] * (double) (n - i + 1) + aa[1] * (double) (n - i) + aa[2] * (double) (n - i - 1);
  265. s = aa[0] + aa[1];
  266. a_inv[i][i] = s * (double) (n - i) + aa[2] * (double) (n - i - 1);
  267. s += aa[2];
  268. for (j = i + 1; j < n; j++) a_inv[i][j] = s * (double) (n - j);
  269. }
  270. //последняя (n-1)
  271. aa[0] = a_inv[n - 1][n - 2];
  272. aa[1] = a_inv[n - 1][n - 1];
  273. s = aa[0] + aa[0] + aa[1];
  274. for (j = 0; j < n - 1; j++) a_inv[n - 1][j] = s;
  275. a_inv[n - 1][n - 1] = aa[0] + aa[1];
  276. break;
  277.  
  278. }//schema
  279. break;
  280.  
  281. case 2: //одна жорданова клетка 2x2 при минимальном с.з.
  282. System.out.println(" J_2 type matrix: must be n > 2");
  283. System.out.println(" schema = " + schema);
  284.  
  285. switch (schema) {
  286. case 1:
  287. //TJ
  288. //первая строка
  289. a[0][0] = J[0];
  290. a[0][1] = 1. - J[0];
  291. //вторая строка
  292. a[1][0] = -J[0];
  293. a[1][1] = -1. + J[0] + J[0];
  294. a[1][2] = -J[2];
  295. //третья строка
  296. a[2][1] = -J[0];
  297. a[2][2] = J[2] + J[2];
  298. if (n > 3) a[2][3] = -J[3];
  299. //до последней
  300. for (i = 3; i < n - 1; i++) {
  301. a[i][i - 1] = -J[i - 1];
  302. a[i][i] = J[i] + J[i];
  303. a[i][i + 1] = -J[i + 1];
  304. }
  305. //последняя (n-1)
  306. if (n > 3) {
  307. a[n - 1][n - 2] = -J[n - 2];
  308. a[n - 1][n - 1] = J[n - 1] + J[n - 1];
  309. }
  310.  
  311. //(TJ)T^{-1}
  312. //первая строка
  313. aa[1] = a[0][0];
  314. aa[2] = a[0][1];
  315. a[0][0] = aa[1] * (double) n + aa[2] * (double) (n - 1);
  316. double s = aa[1] + aa[2];
  317. for (j = 1; j < n; j++) a[0][j] = s * (double) (n - j);
  318. //до последней
  319. for (i = 1; i < n - 1; i++) {
  320. aa[0] = a[i][i - 1];
  321. aa[1] = a[i][i];
  322. aa[2] = a[i][i + 1];
  323. for (j = 0; j < i; j++)
  324. a[i][j] = aa[0] * (double) (n - i + 1) + aa[1] * (double) (n - i) + aa[2] * (double) (n - i - 1);
  325. s = aa[0] + aa[1];
  326. a[i][i] = s * (double) (n - i) + aa[2] * (double) (n - i - 1);
  327. s += aa[2];
  328. for (j = i + 1; j < n; j++) a[i][j] = s * (double) (n - j);
  329. }
  330. //последняя (n-1)
  331. aa[0] = a[n - 1][n - 2];
  332. aa[1] = a[n - 1][n - 1];
  333. s = aa[0] + aa[0] + aa[1];
  334. for (j = 0; j < n - 1; j++) a[n - 1][j] = s;
  335. a[n - 1][n - 1] = aa[0] + aa[1];
  336. //_______
  337. //TJ^{-1}
  338. //первая строка
  339. a_inv[0][0] = J_inv[0];
  340. a_inv[0][1] = -J_inv[0] * J_inv[0] - J_inv[0];
  341. //вторая строка
  342. a_inv[1][0] = -J_inv[0];
  343. a_inv[1][1] = J_inv[0] * J_inv[0] + J_inv[0] + J_inv[0];
  344. a_inv[1][2] = -J_inv[2];
  345. //третья строка
  346. a_inv[2][1] = -J_inv[0];
  347. a_inv[2][2] = J_inv[2] + J_inv[2];
  348. if (n > 3) a_inv[2][3] = -J_inv[3];
  349. //до последней
  350. for (i = 3; i < n - 1; i++) {
  351. a_inv[i][i - 1] = -J_inv[i - 1];
  352. a_inv[i][i] = J_inv[i] + J_inv[i];
  353. a_inv[i][i + 1] = -J_inv[i + 1];
  354. }
  355. //последняя (n-1)
  356. if (n > 3) {
  357. a_inv[n - 1][n - 2] = -J_inv[n - 2];
  358. a_inv[n - 1][n - 1] = J_inv[n - 1] + J_inv[n - 1];
  359. }
  360.  
  361. //(TJ^{-1})T^{-1}
  362. //первая строка
  363. aa[1] = a_inv[0][0];
  364. aa[2] = a_inv[0][1];
  365. a_inv[0][0] = aa[1] * (double) n + aa[2] * (double) (n - 1);
  366. s = aa[1] + aa[2];
  367. for (j = 1; j < n; j++) a_inv[0][j] = s * (double) (n - j);
  368. //до последней
  369. for (i = 1; i < n - 1; i++) {
  370. aa[0] = a_inv[i][i - 1];
  371. aa[1] = a_inv[i][i];
  372. aa[2] = a_inv[i][i + 1];
  373. for (j = 0; j < i; j++)
  374. a_inv[i][j] = aa[0] * (double) (n - i + 1) + aa[1] * (double) (n - i) + aa[2] * (double) (n - i - 1);
  375. s = aa[0] + aa[1];
  376. a_inv[i][i] = s * (double) (n - i) + aa[2] * (double) (n - i - 1);
  377. s += aa[2];
  378. for (j = i + 1; j < n; j++) a_inv[i][j] = s * (double) (n - j);
  379. }
  380. //последняя (n-1)
  381. aa[0] = a_inv[n - 1][n - 2];
  382. aa[1] = a_inv[n - 1][n - 1];
  383. s = aa[0] + aa[0] + aa[1];
  384. for (j = 0; j < n - 1; j++) a_inv[n - 1][j] = s;
  385. a_inv[n - 1][n - 1] = aa[0] + aa[1];
  386.  
  387.  
  388. break;
  389. }//schema
  390.  
  391. break;
  392.  
  393. }//variant
  394.  
  395. //______________________________________________________________________
  396.  
  397. double norm, norm_inv;
  398.  
  399. norm = matr_inf_norm(a, n);
  400. System.out.println(" || A || = " + norm);
  401.  
  402. norm_inv = matr_inf_norm(a_inv, n);
  403. System.out.println(" ||A_inv|| = " + norm_inv);
  404.  
  405. double obusl = norm * norm_inv;
  406. System.out.println(" obusl = " + obusl);
  407.  
  408. //невязка генерации
  409. double[][] r = new double[n][];
  410. for (i = 0; i < n; i++)
  411. r[i] = new double[n];
  412. matr_mul(a, a_inv, r, n);
  413. for (i = 0; i < n; i++) r[i][i] -= 1.;
  414.  
  415. /* cout << "r:" << endl;
  416. for( i = 0; i < n; i++ )
  417. {
  418. for( j = 0; j < n; j++ ) cout << " " << r[i][j];
  419. cout << endl;
  420. }
  421. */
  422. norm = matr_inf_norm(r, n);
  423. System.out.println(" ||R_gen|| = " + norm);
  424.  
  425.  
  426. }//mygen
  427.  
  428.  
  429. }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!