[2806] | 1 | /*************************************************************************
|
---|
| 2 | Copyright (c) 2005-2007, Sergey Bochkanov (ALGLIB project).
|
---|
| 3 |
|
---|
| 4 | >>> SOURCE LICENSE >>>
|
---|
| 5 | This program is free software; you can redistribute it and/or modify
|
---|
| 6 | it under the terms of the GNU General Public License as published by
|
---|
| 7 | the Free Software Foundation (www.fsf.org); either version 2 of the
|
---|
| 8 | License, or (at your option) any later version.
|
---|
| 9 |
|
---|
| 10 | This program is distributed in the hope that it will be useful,
|
---|
| 11 | but WITHOUT ANY WARRANTY; without even the implied warranty of
|
---|
| 12 | MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
---|
| 13 | GNU General Public License for more details.
|
---|
| 14 |
|
---|
| 15 | A copy of the GNU General Public License is available at
|
---|
| 16 | http://www.fsf.org/licensing/licenses
|
---|
| 17 |
|
---|
| 18 | >>> END OF LICENSE >>>
|
---|
| 19 | *************************************************************************/
|
---|
| 20 |
|
---|
| 21 | using System;
|
---|
| 22 |
|
---|
| 23 | namespace alglib
|
---|
| 24 | {
|
---|
| 25 | public class hbisinv
|
---|
| 26 | {
|
---|
| 27 | /*************************************************************************
|
---|
| 28 | Subroutine for finding the eigenvalues (and eigenvectors) of a Hermitian
|
---|
| 29 | matrix in a given half-interval (A, B] by using a bisection and inverse
|
---|
| 30 | iteration
|
---|
| 31 |
|
---|
| 32 | Input parameters:
|
---|
| 33 | A - Hermitian matrix which is given by its upper or lower
|
---|
| 34 | triangular part. Array whose indexes range within
|
---|
| 35 | [0..N-1, 0..N-1].
|
---|
| 36 | N - size of matrix A.
|
---|
| 37 | ZNeeded - flag controlling whether the eigenvectors are needed or
|
---|
| 38 | not. If ZNeeded is equal to:
|
---|
| 39 | * 0, the eigenvectors are not returned;
|
---|
| 40 | * 1, the eigenvectors are returned.
|
---|
| 41 | IsUpperA - storage format of matrix A.
|
---|
| 42 | B1, B2 - half-interval (B1, B2] to search eigenvalues in.
|
---|
| 43 |
|
---|
| 44 | Output parameters:
|
---|
| 45 | M - number of eigenvalues found in a given half-interval, M>=0
|
---|
| 46 | W - array of the eigenvalues found.
|
---|
| 47 | Array whose index ranges within [0..M-1].
|
---|
| 48 | Z - if ZNeeded is equal to:
|
---|
| 49 | * 0, Z hasnt changed;
|
---|
| 50 | * 1, Z contains eigenvectors.
|
---|
| 51 | Array whose indexes range within [0..N-1, 0..M-1].
|
---|
| 52 | The eigenvectors are stored in the matrix columns.
|
---|
| 53 |
|
---|
| 54 | Result:
|
---|
| 55 | True, if successful. M contains the number of eigenvalues in the given
|
---|
| 56 | half-interval (could be equal to 0), W contains the eigenvalues,
|
---|
| 57 | Z contains the eigenvectors (if needed).
|
---|
| 58 |
|
---|
| 59 | False, if the bisection method subroutine wasn't able to find the
|
---|
| 60 | eigenvalues in the given interval or if the inverse iteration
|
---|
| 61 | subroutine wasn't able to find all the corresponding eigenvectors.
|
---|
| 62 | In that case, the eigenvalues and eigenvectors are not returned, M is
|
---|
| 63 | equal to 0.
|
---|
| 64 |
|
---|
| 65 | Note:
|
---|
| 66 | eigen vectors of Hermitian matrix are defined up to multiplication by
|
---|
| 67 | a complex number L, such as |L|=1.
|
---|
| 68 |
|
---|
| 69 | -- ALGLIB --
|
---|
| 70 | Copyright 07.01.2006, 24.03.2007 by Bochkanov Sergey.
|
---|
| 71 | *************************************************************************/
|
---|
| 72 | public static bool hmatrixevdr(AP.Complex[,] a,
|
---|
| 73 | int n,
|
---|
| 74 | int zneeded,
|
---|
| 75 | bool isupper,
|
---|
| 76 | double b1,
|
---|
| 77 | double b2,
|
---|
| 78 | ref int m,
|
---|
| 79 | ref double[] w,
|
---|
| 80 | ref AP.Complex[,] z)
|
---|
| 81 | {
|
---|
| 82 | bool result = new bool();
|
---|
| 83 | AP.Complex[,] q = new AP.Complex[0,0];
|
---|
| 84 | double[,] t = new double[0,0];
|
---|
| 85 | AP.Complex[] tau = new AP.Complex[0];
|
---|
| 86 | double[] e = new double[0];
|
---|
| 87 | double[] work = new double[0];
|
---|
| 88 | int i = 0;
|
---|
| 89 | int k = 0;
|
---|
| 90 | double v = 0;
|
---|
| 91 | int i_ = 0;
|
---|
| 92 |
|
---|
| 93 | a = (AP.Complex[,])a.Clone();
|
---|
| 94 |
|
---|
| 95 | System.Diagnostics.Debug.Assert(zneeded==0 | zneeded==1, "HermitianEigenValuesAndVectorsInInterval: incorrect ZNeeded");
|
---|
| 96 |
|
---|
| 97 | //
|
---|
| 98 | // Reduce to tridiagonal form
|
---|
| 99 | //
|
---|
| 100 | htridiagonal.hmatrixtd(ref a, n, isupper, ref tau, ref w, ref e);
|
---|
| 101 | if( zneeded==1 )
|
---|
| 102 | {
|
---|
| 103 | htridiagonal.hmatrixtdunpackq(ref a, n, isupper, ref tau, ref q);
|
---|
| 104 | zneeded = 2;
|
---|
| 105 | }
|
---|
| 106 |
|
---|
| 107 | //
|
---|
| 108 | // Bisection and inverse iteration
|
---|
| 109 | //
|
---|
| 110 | result = tdbisinv.smatrixtdevdr(ref w, ref e, n, zneeded, b1, b2, ref m, ref t);
|
---|
| 111 |
|
---|
| 112 | //
|
---|
| 113 | // Eigenvectors are needed
|
---|
| 114 | // Calculate Z = Q*T = Re(Q)*T + i*Im(Q)*T
|
---|
| 115 | //
|
---|
| 116 | if( result & zneeded!=0 & m!=0 )
|
---|
| 117 | {
|
---|
| 118 | work = new double[m-1+1];
|
---|
| 119 | z = new AP.Complex[n-1+1, m-1+1];
|
---|
| 120 | for(i=0; i<=n-1; i++)
|
---|
| 121 | {
|
---|
| 122 |
|
---|
| 123 | //
|
---|
| 124 | // Calculate real part
|
---|
| 125 | //
|
---|
| 126 | for(k=0; k<=m-1; k++)
|
---|
| 127 | {
|
---|
| 128 | work[k] = 0;
|
---|
| 129 | }
|
---|
| 130 | for(k=0; k<=n-1; k++)
|
---|
| 131 | {
|
---|
| 132 | v = q[i,k].x;
|
---|
| 133 | for(i_=0; i_<=m-1;i_++)
|
---|
| 134 | {
|
---|
| 135 | work[i_] = work[i_] + v*t[k,i_];
|
---|
| 136 | }
|
---|
| 137 | }
|
---|
| 138 | for(k=0; k<=m-1; k++)
|
---|
| 139 | {
|
---|
| 140 | z[i,k].x = work[k];
|
---|
| 141 | }
|
---|
| 142 |
|
---|
| 143 | //
|
---|
| 144 | // Calculate imaginary part
|
---|
| 145 | //
|
---|
| 146 | for(k=0; k<=m-1; k++)
|
---|
| 147 | {
|
---|
| 148 | work[k] = 0;
|
---|
| 149 | }
|
---|
| 150 | for(k=0; k<=n-1; k++)
|
---|
| 151 | {
|
---|
| 152 | v = q[i,k].y;
|
---|
| 153 | for(i_=0; i_<=m-1;i_++)
|
---|
| 154 | {
|
---|
| 155 | work[i_] = work[i_] + v*t[k,i_];
|
---|
| 156 | }
|
---|
| 157 | }
|
---|
| 158 | for(k=0; k<=m-1; k++)
|
---|
| 159 | {
|
---|
| 160 | z[i,k].y = work[k];
|
---|
| 161 | }
|
---|
| 162 | }
|
---|
| 163 | }
|
---|
| 164 | return result;
|
---|
| 165 | }
|
---|
| 166 |
|
---|
| 167 |
|
---|
| 168 | /*************************************************************************
|
---|
| 169 | Subroutine for finding the eigenvalues and eigenvectors of a Hermitian
|
---|
| 170 | matrix with given indexes by using bisection and inverse iteration methods
|
---|
| 171 |
|
---|
| 172 | Input parameters:
|
---|
| 173 | A - Hermitian matrix which is given by its upper or lower
|
---|
| 174 | triangular part.
|
---|
| 175 | Array whose indexes range within [0..N-1, 0..N-1].
|
---|
| 176 | N - size of matrix A.
|
---|
| 177 | ZNeeded - flag controlling whether the eigenvectors are needed or
|
---|
| 178 | not. If ZNeeded is equal to:
|
---|
| 179 | * 0, the eigenvectors are not returned;
|
---|
| 180 | * 1, the eigenvectors are returned.
|
---|
| 181 | IsUpperA - storage format of matrix A.
|
---|
| 182 | I1, I2 - index interval for searching (from I1 to I2).
|
---|
| 183 | 0 <= I1 <= I2 <= N-1.
|
---|
| 184 |
|
---|
| 185 | Output parameters:
|
---|
| 186 | W - array of the eigenvalues found.
|
---|
| 187 | Array whose index ranges within [0..I2-I1].
|
---|
| 188 | Z - if ZNeeded is equal to:
|
---|
| 189 | * 0, Z hasnt changed;
|
---|
| 190 | * 1, Z contains eigenvectors.
|
---|
| 191 | Array whose indexes range within [0..N-1, 0..I2-I1].
|
---|
| 192 | In that case, the eigenvectors are stored in the matrix
|
---|
| 193 | columns.
|
---|
| 194 |
|
---|
| 195 | Result:
|
---|
| 196 | True, if successful. W contains the eigenvalues, Z contains the
|
---|
| 197 | eigenvectors (if needed).
|
---|
| 198 |
|
---|
| 199 | False, if the bisection method subroutine wasn't able to find the
|
---|
| 200 | eigenvalues in the given interval or if the inverse iteration
|
---|
| 201 | subroutine wasn't able to find all the corresponding eigenvectors.
|
---|
| 202 | In that case, the eigenvalues and eigenvectors are not returned.
|
---|
| 203 |
|
---|
| 204 | Note:
|
---|
| 205 | eigen vectors of Hermitian matrix are defined up to multiplication by
|
---|
| 206 | a complex number L, such as |L|=1.
|
---|
| 207 |
|
---|
| 208 | -- ALGLIB --
|
---|
| 209 | Copyright 07.01.2006, 24.03.2007 by Bochkanov Sergey.
|
---|
| 210 | *************************************************************************/
|
---|
| 211 | public static bool hmatrixevdi(AP.Complex[,] a,
|
---|
| 212 | int n,
|
---|
| 213 | int zneeded,
|
---|
| 214 | bool isupper,
|
---|
| 215 | int i1,
|
---|
| 216 | int i2,
|
---|
| 217 | ref double[] w,
|
---|
| 218 | ref AP.Complex[,] z)
|
---|
| 219 | {
|
---|
| 220 | bool result = new bool();
|
---|
| 221 | AP.Complex[,] q = new AP.Complex[0,0];
|
---|
| 222 | double[,] t = new double[0,0];
|
---|
| 223 | AP.Complex[] tau = new AP.Complex[0];
|
---|
| 224 | double[] e = new double[0];
|
---|
| 225 | double[] work = new double[0];
|
---|
| 226 | int i = 0;
|
---|
| 227 | int k = 0;
|
---|
| 228 | double v = 0;
|
---|
| 229 | int m = 0;
|
---|
| 230 | int i_ = 0;
|
---|
| 231 |
|
---|
| 232 | a = (AP.Complex[,])a.Clone();
|
---|
| 233 |
|
---|
| 234 | System.Diagnostics.Debug.Assert(zneeded==0 | zneeded==1, "HermitianEigenValuesAndVectorsByIndexes: incorrect ZNeeded");
|
---|
| 235 |
|
---|
| 236 | //
|
---|
| 237 | // Reduce to tridiagonal form
|
---|
| 238 | //
|
---|
| 239 | htridiagonal.hmatrixtd(ref a, n, isupper, ref tau, ref w, ref e);
|
---|
| 240 | if( zneeded==1 )
|
---|
| 241 | {
|
---|
| 242 | htridiagonal.hmatrixtdunpackq(ref a, n, isupper, ref tau, ref q);
|
---|
| 243 | zneeded = 2;
|
---|
| 244 | }
|
---|
| 245 |
|
---|
| 246 | //
|
---|
| 247 | // Bisection and inverse iteration
|
---|
| 248 | //
|
---|
| 249 | result = tdbisinv.smatrixtdevdi(ref w, ref e, n, zneeded, i1, i2, ref t);
|
---|
| 250 |
|
---|
| 251 | //
|
---|
| 252 | // Eigenvectors are needed
|
---|
| 253 | // Calculate Z = Q*T = Re(Q)*T + i*Im(Q)*T
|
---|
| 254 | //
|
---|
| 255 | m = i2-i1+1;
|
---|
| 256 | if( result & zneeded!=0 )
|
---|
| 257 | {
|
---|
| 258 | work = new double[m-1+1];
|
---|
| 259 | z = new AP.Complex[n-1+1, m-1+1];
|
---|
| 260 | for(i=0; i<=n-1; i++)
|
---|
| 261 | {
|
---|
| 262 |
|
---|
| 263 | //
|
---|
| 264 | // Calculate real part
|
---|
| 265 | //
|
---|
| 266 | for(k=0; k<=m-1; k++)
|
---|
| 267 | {
|
---|
| 268 | work[k] = 0;
|
---|
| 269 | }
|
---|
| 270 | for(k=0; k<=n-1; k++)
|
---|
| 271 | {
|
---|
| 272 | v = q[i,k].x;
|
---|
| 273 | for(i_=0; i_<=m-1;i_++)
|
---|
| 274 | {
|
---|
| 275 | work[i_] = work[i_] + v*t[k,i_];
|
---|
| 276 | }
|
---|
| 277 | }
|
---|
| 278 | for(k=0; k<=m-1; k++)
|
---|
| 279 | {
|
---|
| 280 | z[i,k].x = work[k];
|
---|
| 281 | }
|
---|
| 282 |
|
---|
| 283 | //
|
---|
| 284 | // Calculate imaginary part
|
---|
| 285 | //
|
---|
| 286 | for(k=0; k<=m-1; k++)
|
---|
| 287 | {
|
---|
| 288 | work[k] = 0;
|
---|
| 289 | }
|
---|
| 290 | for(k=0; k<=n-1; k++)
|
---|
| 291 | {
|
---|
| 292 | v = q[i,k].y;
|
---|
| 293 | for(i_=0; i_<=m-1;i_++)
|
---|
| 294 | {
|
---|
| 295 | work[i_] = work[i_] + v*t[k,i_];
|
---|
| 296 | }
|
---|
| 297 | }
|
---|
| 298 | for(k=0; k<=m-1; k++)
|
---|
| 299 | {
|
---|
| 300 | z[i,k].y = work[k];
|
---|
| 301 | }
|
---|
| 302 | }
|
---|
| 303 | }
|
---|
| 304 | return result;
|
---|
| 305 | }
|
---|
| 306 |
|
---|
| 307 |
|
---|
| 308 | public static bool hermitianeigenvaluesandvectorsininterval(AP.Complex[,] a,
|
---|
| 309 | int n,
|
---|
| 310 | int zneeded,
|
---|
| 311 | bool isupper,
|
---|
| 312 | double b1,
|
---|
| 313 | double b2,
|
---|
| 314 | ref int m,
|
---|
| 315 | ref double[] w,
|
---|
| 316 | ref AP.Complex[,] z)
|
---|
| 317 | {
|
---|
| 318 | bool result = new bool();
|
---|
| 319 | AP.Complex[,] q = new AP.Complex[0,0];
|
---|
| 320 | double[,] t = new double[0,0];
|
---|
| 321 | AP.Complex[] tau = new AP.Complex[0];
|
---|
| 322 | double[] e = new double[0];
|
---|
| 323 | double[] work = new double[0];
|
---|
| 324 | int i = 0;
|
---|
| 325 | int k = 0;
|
---|
| 326 | double v = 0;
|
---|
| 327 | int i_ = 0;
|
---|
| 328 |
|
---|
| 329 | a = (AP.Complex[,])a.Clone();
|
---|
| 330 |
|
---|
| 331 | System.Diagnostics.Debug.Assert(zneeded==0 | zneeded==1, "HermitianEigenValuesAndVectorsInInterval: incorrect ZNeeded");
|
---|
| 332 |
|
---|
| 333 | //
|
---|
| 334 | // Reduce to tridiagonal form
|
---|
| 335 | //
|
---|
| 336 | htridiagonal.hermitiantotridiagonal(ref a, n, isupper, ref tau, ref w, ref e);
|
---|
| 337 | if( zneeded==1 )
|
---|
| 338 | {
|
---|
| 339 | htridiagonal.unpackqfromhermitiantridiagonal(ref a, n, isupper, ref tau, ref q);
|
---|
| 340 | zneeded = 2;
|
---|
| 341 | }
|
---|
| 342 |
|
---|
| 343 | //
|
---|
| 344 | // Bisection and inverse iteration
|
---|
| 345 | //
|
---|
| 346 | result = tdbisinv.tridiagonaleigenvaluesandvectorsininterval(ref w, ref e, n, zneeded, b1, b2, ref m, ref t);
|
---|
| 347 |
|
---|
| 348 | //
|
---|
| 349 | // Eigenvectors are needed
|
---|
| 350 | // Calculate Z = Q*T = Re(Q)*T + i*Im(Q)*T
|
---|
| 351 | //
|
---|
| 352 | if( result & zneeded!=0 & m!=0 )
|
---|
| 353 | {
|
---|
| 354 | work = new double[m+1];
|
---|
| 355 | z = new AP.Complex[n+1, m+1];
|
---|
| 356 | for(i=1; i<=n; i++)
|
---|
| 357 | {
|
---|
| 358 |
|
---|
| 359 | //
|
---|
| 360 | // Calculate real part
|
---|
| 361 | //
|
---|
| 362 | for(k=1; k<=m; k++)
|
---|
| 363 | {
|
---|
| 364 | work[k] = 0;
|
---|
| 365 | }
|
---|
| 366 | for(k=1; k<=n; k++)
|
---|
| 367 | {
|
---|
| 368 | v = q[i,k].x;
|
---|
| 369 | for(i_=1; i_<=m;i_++)
|
---|
| 370 | {
|
---|
| 371 | work[i_] = work[i_] + v*t[k,i_];
|
---|
| 372 | }
|
---|
| 373 | }
|
---|
| 374 | for(k=1; k<=m; k++)
|
---|
| 375 | {
|
---|
| 376 | z[i,k].x = work[k];
|
---|
| 377 | }
|
---|
| 378 |
|
---|
| 379 | //
|
---|
| 380 | // Calculate imaginary part
|
---|
| 381 | //
|
---|
| 382 | for(k=1; k<=m; k++)
|
---|
| 383 | {
|
---|
| 384 | work[k] = 0;
|
---|
| 385 | }
|
---|
| 386 | for(k=1; k<=n; k++)
|
---|
| 387 | {
|
---|
| 388 | v = q[i,k].y;
|
---|
| 389 | for(i_=1; i_<=m;i_++)
|
---|
| 390 | {
|
---|
| 391 | work[i_] = work[i_] + v*t[k,i_];
|
---|
| 392 | }
|
---|
| 393 | }
|
---|
| 394 | for(k=1; k<=m; k++)
|
---|
| 395 | {
|
---|
| 396 | z[i,k].y = work[k];
|
---|
| 397 | }
|
---|
| 398 | }
|
---|
| 399 | }
|
---|
| 400 | return result;
|
---|
| 401 | }
|
---|
| 402 |
|
---|
| 403 |
|
---|
| 404 | public static bool hermitianeigenvaluesandvectorsbyindexes(AP.Complex[,] a,
|
---|
| 405 | int n,
|
---|
| 406 | int zneeded,
|
---|
| 407 | bool isupper,
|
---|
| 408 | int i1,
|
---|
| 409 | int i2,
|
---|
| 410 | ref double[] w,
|
---|
| 411 | ref AP.Complex[,] z)
|
---|
| 412 | {
|
---|
| 413 | bool result = new bool();
|
---|
| 414 | AP.Complex[,] q = new AP.Complex[0,0];
|
---|
| 415 | double[,] t = new double[0,0];
|
---|
| 416 | AP.Complex[] tau = new AP.Complex[0];
|
---|
| 417 | double[] e = new double[0];
|
---|
| 418 | double[] work = new double[0];
|
---|
| 419 | int i = 0;
|
---|
| 420 | int k = 0;
|
---|
| 421 | double v = 0;
|
---|
| 422 | int m = 0;
|
---|
| 423 | int i_ = 0;
|
---|
| 424 |
|
---|
| 425 | a = (AP.Complex[,])a.Clone();
|
---|
| 426 |
|
---|
| 427 | System.Diagnostics.Debug.Assert(zneeded==0 | zneeded==1, "HermitianEigenValuesAndVectorsByIndexes: incorrect ZNeeded");
|
---|
| 428 |
|
---|
| 429 | //
|
---|
| 430 | // Reduce to tridiagonal form
|
---|
| 431 | //
|
---|
| 432 | htridiagonal.hermitiantotridiagonal(ref a, n, isupper, ref tau, ref w, ref e);
|
---|
| 433 | if( zneeded==1 )
|
---|
| 434 | {
|
---|
| 435 | htridiagonal.unpackqfromhermitiantridiagonal(ref a, n, isupper, ref tau, ref q);
|
---|
| 436 | zneeded = 2;
|
---|
| 437 | }
|
---|
| 438 |
|
---|
| 439 | //
|
---|
| 440 | // Bisection and inverse iteration
|
---|
| 441 | //
|
---|
| 442 | result = tdbisinv.tridiagonaleigenvaluesandvectorsbyindexes(ref w, ref e, n, zneeded, i1, i2, ref t);
|
---|
| 443 |
|
---|
| 444 | //
|
---|
| 445 | // Eigenvectors are needed
|
---|
| 446 | // Calculate Z = Q*T = Re(Q)*T + i*Im(Q)*T
|
---|
| 447 | //
|
---|
| 448 | m = i2-i1+1;
|
---|
| 449 | if( result & zneeded!=0 )
|
---|
| 450 | {
|
---|
| 451 | work = new double[m+1];
|
---|
| 452 | z = new AP.Complex[n+1, m+1];
|
---|
| 453 | for(i=1; i<=n; i++)
|
---|
| 454 | {
|
---|
| 455 |
|
---|
| 456 | //
|
---|
| 457 | // Calculate real part
|
---|
| 458 | //
|
---|
| 459 | for(k=1; k<=m; k++)
|
---|
| 460 | {
|
---|
| 461 | work[k] = 0;
|
---|
| 462 | }
|
---|
| 463 | for(k=1; k<=n; k++)
|
---|
| 464 | {
|
---|
| 465 | v = q[i,k].x;
|
---|
| 466 | for(i_=1; i_<=m;i_++)
|
---|
| 467 | {
|
---|
| 468 | work[i_] = work[i_] + v*t[k,i_];
|
---|
| 469 | }
|
---|
| 470 | }
|
---|
| 471 | for(k=1; k<=m; k++)
|
---|
| 472 | {
|
---|
| 473 | z[i,k].x = work[k];
|
---|
| 474 | }
|
---|
| 475 |
|
---|
| 476 | //
|
---|
| 477 | // Calculate imaginary part
|
---|
| 478 | //
|
---|
| 479 | for(k=1; k<=m; k++)
|
---|
| 480 | {
|
---|
| 481 | work[k] = 0;
|
---|
| 482 | }
|
---|
| 483 | for(k=1; k<=n; k++)
|
---|
| 484 | {
|
---|
| 485 | v = q[i,k].y;
|
---|
| 486 | for(i_=1; i_<=m;i_++)
|
---|
| 487 | {
|
---|
| 488 | work[i_] = work[i_] + v*t[k,i_];
|
---|
| 489 | }
|
---|
| 490 | }
|
---|
| 491 | for(k=1; k<=m; k++)
|
---|
| 492 | {
|
---|
| 493 | z[i,k].y = work[k];
|
---|
| 494 | }
|
---|
| 495 | }
|
---|
| 496 | }
|
---|
| 497 | return result;
|
---|
| 498 | }
|
---|
| 499 | }
|
---|
| 500 | }
|
---|