1 2//--------------------------------------------------------------------- 3//--------------------------------------------------------------------- 4#define ablock(m,n) ablock[(m-1)+5*(n-1)] 5#define bblock(m,n) bblock[(m-1)+5*(n-1)] 6#define cblock(m,n) cblock[(m-1)+5*(n-1)] 7#define avec(m) avec[m-1] 8#define bvec(m) bvec[m-1] 9#define lhs(m,n) lhs[(m-1)+5*(n-1)] 10#define c(m,n) c[(m-1)+5*(n-1)] 11#define r(m) r[m-1] 12 13void matvec_sub(double ablock[],double avec[],double bvec[]) { 14 15//--------------------------------------------------------------------- 16//--------------------------------------------------------------------- 17 18//--------------------------------------------------------------------- 19// subtracts bvec=bvec - ablock*avec 20//--------------------------------------------------------------------- 21 22//--------------------------------------------------------------------- 23// rhs(i,ic,jc,kc,ccell) = rhs(i,ic,jc,kc,ccell) 24// $ - lhs(i,1,ablock,ia,ja,ka,acell)* 25//--------------------------------------------------------------------- 26 bvec(1) = bvec(1) - ablock(1,1)*avec(1) 27 - ablock(1,2)*avec(2) 28 - ablock(1,3)*avec(3) 29 - ablock(1,4)*avec(4) 30 - ablock(1,5)*avec(5); 31 bvec(2) = bvec(2) - ablock(2,1)*avec(1) 32 - ablock(2,2)*avec(2) 33 - ablock(2,3)*avec(3) 34 - ablock(2,4)*avec(4) 35 - ablock(2,5)*avec(5); 36 bvec(3) = bvec(3) - ablock(3,1)*avec(1) 37 - ablock(3,2)*avec(2) 38 - ablock(3,3)*avec(3) 39 - ablock(3,4)*avec(4) 40 - ablock(3,5)*avec(5); 41 bvec(4) = bvec(4) - ablock(4,1)*avec(1) 42 - ablock(4,2)*avec(2) 43 - ablock(4,3)*avec(3) 44 - ablock(4,4)*avec(4) 45 - ablock(4,5)*avec(5); 46 bvec(5) = bvec(5) - ablock(5,1)*avec(1) 47 - ablock(5,2)*avec(2) 48 - ablock(5,3)*avec(3) 49 - ablock(5,4)*avec(4) 50 - ablock(5,5)*avec(5); 51 52 53 return; 54} 55 56//--------------------------------------------------------------------- 57//--------------------------------------------------------------------- 58 59void matmul_sub(double ablock[], double bblock[], double cblock[]) { 60 61//--------------------------------------------------------------------- 62//--------------------------------------------------------------------- 63 64//--------------------------------------------------------------------- 65// subtracts a(i,j,k) X b(i,j,k) from c(i,j,k) 66//--------------------------------------------------------------------- 67 68 69 cblock(1,1) = cblock(1,1) - ablock(1,1)*bblock(1,1) 70 - ablock(1,2)*bblock(2,1) 71 - ablock(1,3)*bblock(3,1) 72 - ablock(1,4)*bblock(4,1) 73 - ablock(1,5)*bblock(5,1); 74 cblock(2,1) = cblock(2,1) - ablock(2,1)*bblock(1,1) 75 - ablock(2,2)*bblock(2,1) 76 - ablock(2,3)*bblock(3,1) 77 - ablock(2,4)*bblock(4,1) 78 - ablock(2,5)*bblock(5,1); 79 cblock(3,1) = cblock(3,1) - ablock(3,1)*bblock(1,1) 80 - ablock(3,2)*bblock(2,1) 81 - ablock(3,3)*bblock(3,1) 82 - ablock(3,4)*bblock(4,1) 83 - ablock(3,5)*bblock(5,1); 84 cblock(4,1) = cblock(4,1) - ablock(4,1)*bblock(1,1) 85 - ablock(4,2)*bblock(2,1) 86 - ablock(4,3)*bblock(3,1) 87 - ablock(4,4)*bblock(4,1) 88 - ablock(4,5)*bblock(5,1); 89 cblock(5,1) = cblock(5,1) - ablock(5,1)*bblock(1,1) 90 - ablock(5,2)*bblock(2,1) 91 - ablock(5,3)*bblock(3,1) 92 - ablock(5,4)*bblock(4,1) 93 - ablock(5,5)*bblock(5,1); 94 cblock(1,2) = cblock(1,2) - ablock(1,1)*bblock(1,2) 95 - ablock(1,2)*bblock(2,2) 96 - ablock(1,3)*bblock(3,2) 97 - ablock(1,4)*bblock(4,2) 98 - ablock(1,5)*bblock(5,2); 99 cblock(2,2) = cblock(2,2) - ablock(2,1)*bblock(1,2) 100 - ablock(2,2)*bblock(2,2) 101 - ablock(2,3)*bblock(3,2) 102 - ablock(2,4)*bblock(4,2) 103 - ablock(2,5)*bblock(5,2); 104 cblock(3,2) = cblock(3,2) - ablock(3,1)*bblock(1,2) 105 - ablock(3,2)*bblock(2,2) 106 - ablock(3,3)*bblock(3,2) 107 - ablock(3,4)*bblock(4,2) 108 - ablock(3,5)*bblock(5,2); 109 cblock(4,2) = cblock(4,2) - ablock(4,1)*bblock(1,2) 110 - ablock(4,2)*bblock(2,2) 111 - ablock(4,3)*bblock(3,2) 112 - ablock(4,4)*bblock(4,2) 113 - ablock(4,5)*bblock(5,2); 114 cblock(5,2) = cblock(5,2) - ablock(5,1)*bblock(1,2) 115 - ablock(5,2)*bblock(2,2) 116 - ablock(5,3)*bblock(3,2) 117 - ablock(5,4)*bblock(4,2) 118 - ablock(5,5)*bblock(5,2); 119 cblock(1,3) = cblock(1,3) - ablock(1,1)*bblock(1,3) 120 - ablock(1,2)*bblock(2,3) 121 - ablock(1,3)*bblock(3,3) 122 - ablock(1,4)*bblock(4,3) 123 - ablock(1,5)*bblock(5,3); 124 cblock(2,3) = cblock(2,3) - ablock(2,1)*bblock(1,3) 125 - ablock(2,2)*bblock(2,3) 126 - ablock(2,3)*bblock(3,3) 127 - ablock(2,4)*bblock(4,3) 128 - ablock(2,5)*bblock(5,3); 129 cblock(3,3) = cblock(3,3) - ablock(3,1)*bblock(1,3) 130 - ablock(3,2)*bblock(2,3) 131 - ablock(3,3)*bblock(3,3) 132 - ablock(3,4)*bblock(4,3) 133 - ablock(3,5)*bblock(5,3); 134 cblock(4,3) = cblock(4,3) - ablock(4,1)*bblock(1,3) 135 - ablock(4,2)*bblock(2,3) 136 - ablock(4,3)*bblock(3,3) 137 - ablock(4,4)*bblock(4,3) 138 - ablock(4,5)*bblock(5,3); 139 cblock(5,3) = cblock(5,3) - ablock(5,1)*bblock(1,3) 140 - ablock(5,2)*bblock(2,3) 141 - ablock(5,3)*bblock(3,3) 142 - ablock(5,4)*bblock(4,3) 143 - ablock(5,5)*bblock(5,3); 144 cblock(1,4) = cblock(1,4) - ablock(1,1)*bblock(1,4) 145 - ablock(1,2)*bblock(2,4) 146 - ablock(1,3)*bblock(3,4) 147 - ablock(1,4)*bblock(4,4) 148 - ablock(1,5)*bblock(5,4); 149 cblock(2,4) = cblock(2,4) - ablock(2,1)*bblock(1,4) 150 - ablock(2,2)*bblock(2,4) 151 - ablock(2,3)*bblock(3,4) 152 - ablock(2,4)*bblock(4,4) 153 - ablock(2,5)*bblock(5,4); 154 cblock(3,4) = cblock(3,4) - ablock(3,1)*bblock(1,4) 155 - ablock(3,2)*bblock(2,4) 156 - ablock(3,3)*bblock(3,4) 157 - ablock(3,4)*bblock(4,4) 158 - ablock(3,5)*bblock(5,4); 159 cblock(4,4) = cblock(4,4) - ablock(4,1)*bblock(1,4) 160 - ablock(4,2)*bblock(2,4) 161 - ablock(4,3)*bblock(3,4) 162 - ablock(4,4)*bblock(4,4) 163 - ablock(4,5)*bblock(5,4); 164 cblock(5,4) = cblock(5,4) - ablock(5,1)*bblock(1,4) 165 - ablock(5,2)*bblock(2,4) 166 - ablock(5,3)*bblock(3,4) 167 - ablock(5,4)*bblock(4,4) 168 - ablock(5,5)*bblock(5,4); 169 cblock(1,5) = cblock(1,5) - ablock(1,1)*bblock(1,5) 170 - ablock(1,2)*bblock(2,5) 171 - ablock(1,3)*bblock(3,5) 172 - ablock(1,4)*bblock(4,5) 173 - ablock(1,5)*bblock(5,5); 174 cblock(2,5) = cblock(2,5) - ablock(2,1)*bblock(1,5) 175 - ablock(2,2)*bblock(2,5) 176 - ablock(2,3)*bblock(3,5) 177 - ablock(2,4)*bblock(4,5) 178 - ablock(2,5)*bblock(5,5); 179 cblock(3,5) = cblock(3,5) - ablock(3,1)*bblock(1,5) 180 - ablock(3,2)*bblock(2,5) 181 - ablock(3,3)*bblock(3,5) 182 - ablock(3,4)*bblock(4,5) 183 - ablock(3,5)*bblock(5,5); 184 cblock(4,5) = cblock(4,5) - ablock(4,1)*bblock(1,5) 185 - ablock(4,2)*bblock(2,5) 186 - ablock(4,3)*bblock(3,5) 187 - ablock(4,4)*bblock(4,5) 188 - ablock(4,5)*bblock(5,5); 189 cblock(5,5) = cblock(5,5) - ablock(5,1)*bblock(1,5) 190 - ablock(5,2)*bblock(2,5) 191 - ablock(5,3)*bblock(3,5) 192 - ablock(5,4)*bblock(4,5) 193 - ablock(5,5)*bblock(5,5); 194 195 196 return; 197} 198 199 200 201//--------------------------------------------------------------------- 202//--------------------------------------------------------------------- 203 204void binvcrhs( double lhs[],double c[],double r[] ) { 205 206//--------------------------------------------------------------------- 207//--------------------------------------------------------------------- 208 209//--------------------------------------------------------------------- 210// 211//--------------------------------------------------------------------- 212 213 double pivot, coeff; 214 215//--------------------------------------------------------------------- 216// 217//--------------------------------------------------------------------- 218 219 pivot = 1.00e0/lhs(1,1); 220 lhs(1,2) = lhs(1,2)*pivot; 221 lhs(1,3) = lhs(1,3)*pivot; 222 lhs(1,4) = lhs(1,4)*pivot; 223 lhs(1,5) = lhs(1,5)*pivot; 224 c(1,1) = c(1,1)*pivot; 225 c(1,2) = c(1,2)*pivot; 226 c(1,3) = c(1,3)*pivot; 227 c(1,4) = c(1,4)*pivot; 228 c(1,5) = c(1,5)*pivot; 229 r(1) = r(1) *pivot; 230 231 coeff = lhs(2,1); 232 lhs(2,2)= lhs(2,2) - coeff*lhs(1,2); 233 lhs(2,3)= lhs(2,3) - coeff*lhs(1,3); 234 lhs(2,4)= lhs(2,4) - coeff*lhs(1,4); 235 lhs(2,5)= lhs(2,5) - coeff*lhs(1,5); 236 c(2,1) = c(2,1) - coeff*c(1,1); 237 c(2,2) = c(2,2) - coeff*c(1,2); 238 c(2,3) = c(2,3) - coeff*c(1,3); 239 c(2,4) = c(2,4) - coeff*c(1,4); 240 c(2,5) = c(2,5) - coeff*c(1,5); 241 r(2) = r(2) - coeff*r(1); 242 243 coeff = lhs(3,1); 244 lhs(3,2)= lhs(3,2) - coeff*lhs(1,2); 245 lhs(3,3)= lhs(3,3) - coeff*lhs(1,3); 246 lhs(3,4)= lhs(3,4) - coeff*lhs(1,4); 247 lhs(3,5)= lhs(3,5) - coeff*lhs(1,5); 248 c(3,1) = c(3,1) - coeff*c(1,1); 249 c(3,2) = c(3,2) - coeff*c(1,2); 250 c(3,3) = c(3,3) - coeff*c(1,3); 251 c(3,4) = c(3,4) - coeff*c(1,4); 252 c(3,5) = c(3,5) - coeff*c(1,5); 253 r(3) = r(3) - coeff*r(1); 254 255 coeff = lhs(4,1); 256 lhs(4,2)= lhs(4,2) - coeff*lhs(1,2); 257 lhs(4,3)= lhs(4,3) - coeff*lhs(1,3); 258 lhs(4,4)= lhs(4,4) - coeff*lhs(1,4); 259 lhs(4,5)= lhs(4,5) - coeff*lhs(1,5); 260 c(4,1) = c(4,1) - coeff*c(1,1); 261 c(4,2) = c(4,2) - coeff*c(1,2); 262 c(4,3) = c(4,3) - coeff*c(1,3); 263 c(4,4) = c(4,4) - coeff*c(1,4); 264 c(4,5) = c(4,5) - coeff*c(1,5); 265 r(4) = r(4) - coeff*r(1); 266 267 coeff = lhs(5,1); 268 lhs(5,2)= lhs(5,2) - coeff*lhs(1,2); 269 lhs(5,3)= lhs(5,3) - coeff*lhs(1,3); 270 lhs(5,4)= lhs(5,4) - coeff*lhs(1,4); 271 lhs(5,5)= lhs(5,5) - coeff*lhs(1,5); 272 c(5,1) = c(5,1) - coeff*c(1,1); 273 c(5,2) = c(5,2) - coeff*c(1,2); 274 c(5,3) = c(5,3) - coeff*c(1,3); 275 c(5,4) = c(5,4) - coeff*c(1,4); 276 c(5,5) = c(5,5) - coeff*c(1,5); 277 r(5) = r(5) - coeff*r(1); 278 279 280 pivot = 1.00e0/lhs(2,2); 281 lhs(2,3) = lhs(2,3)*pivot; 282 lhs(2,4) = lhs(2,4)*pivot; 283 lhs(2,5) = lhs(2,5)*pivot; 284 c(2,1) = c(2,1)*pivot; 285 c(2,2) = c(2,2)*pivot; 286 c(2,3) = c(2,3)*pivot; 287 c(2,4) = c(2,4)*pivot; 288 c(2,5) = c(2,5)*pivot; 289 r(2) = r(2) *pivot; 290 291 coeff = lhs(1,2); 292 lhs(1,3)= lhs(1,3) - coeff*lhs(2,3); 293 lhs(1,4)= lhs(1,4) - coeff*lhs(2,4); 294 lhs(1,5)= lhs(1,5) - coeff*lhs(2,5); 295 c(1,1) = c(1,1) - coeff*c(2,1); 296 c(1,2) = c(1,2) - coeff*c(2,2); 297 c(1,3) = c(1,3) - coeff*c(2,3); 298 c(1,4) = c(1,4) - coeff*c(2,4); 299 c(1,5) = c(1,5) - coeff*c(2,5); 300 r(1) = r(1) - coeff*r(2); 301 302 coeff = lhs(3,2); 303 lhs(3,3)= lhs(3,3) - coeff*lhs(2,3); 304 lhs(3,4)= lhs(3,4) - coeff*lhs(2,4); 305 lhs(3,5)= lhs(3,5) - coeff*lhs(2,5); 306 c(3,1) = c(3,1) - coeff*c(2,1); 307 c(3,2) = c(3,2) - coeff*c(2,2); 308 c(3,3) = c(3,3) - coeff*c(2,3); 309 c(3,4) = c(3,4) - coeff*c(2,4); 310 c(3,5) = c(3,5) - coeff*c(2,5); 311 r(3) = r(3) - coeff*r(2); 312 313 coeff = lhs(4,2); 314 lhs(4,3)= lhs(4,3) - coeff*lhs(2,3); 315 lhs(4,4)= lhs(4,4) - coeff*lhs(2,4); 316 lhs(4,5)= lhs(4,5) - coeff*lhs(2,5); 317 c(4,1) = c(4,1) - coeff*c(2,1); 318 c(4,2) = c(4,2) - coeff*c(2,2); 319 c(4,3) = c(4,3) - coeff*c(2,3); 320 c(4,4) = c(4,4) - coeff*c(2,4); 321 c(4,5) = c(4,5) - coeff*c(2,5); 322 r(4) = r(4) - coeff*r(2); 323 324 coeff = lhs(5,2); 325 lhs(5,3)= lhs(5,3) - coeff*lhs(2,3); 326 lhs(5,4)= lhs(5,4) - coeff*lhs(2,4); 327 lhs(5,5)= lhs(5,5) - coeff*lhs(2,5); 328 c(5,1) = c(5,1) - coeff*c(2,1); 329 c(5,2) = c(5,2) - coeff*c(2,2); 330 c(5,3) = c(5,3) - coeff*c(2,3); 331 c(5,4) = c(5,4) - coeff*c(2,4); 332 c(5,5) = c(5,5) - coeff*c(2,5); 333 r(5) = r(5) - coeff*r(2); 334 335 336 pivot = 1.00e0/lhs(3,3); 337 lhs(3,4) = lhs(3,4)*pivot; 338 lhs(3,5) = lhs(3,5)*pivot; 339 c(3,1) = c(3,1)*pivot; 340 c(3,2) = c(3,2)*pivot; 341 c(3,3) = c(3,3)*pivot; 342 c(3,4) = c(3,4)*pivot; 343 c(3,5) = c(3,5)*pivot; 344 r(3) = r(3) *pivot; 345 346 coeff = lhs(1,3); 347 lhs(1,4)= lhs(1,4) - coeff*lhs(3,4); 348 lhs(1,5)= lhs(1,5) - coeff*lhs(3,5); 349 c(1,1) = c(1,1) - coeff*c(3,1); 350 c(1,2) = c(1,2) - coeff*c(3,2); 351 c(1,3) = c(1,3) - coeff*c(3,3); 352 c(1,4) = c(1,4) - coeff*c(3,4); 353 c(1,5) = c(1,5) - coeff*c(3,5); 354 r(1) = r(1) - coeff*r(3); 355 356 coeff = lhs(2,3); 357 lhs(2,4)= lhs(2,4) - coeff*lhs(3,4); 358 lhs(2,5)= lhs(2,5) - coeff*lhs(3,5); 359 c(2,1) = c(2,1) - coeff*c(3,1); 360 c(2,2) = c(2,2) - coeff*c(3,2); 361 c(2,3) = c(2,3) - coeff*c(3,3); 362 c(2,4) = c(2,4) - coeff*c(3,4); 363 c(2,5) = c(2,5) - coeff*c(3,5); 364 r(2) = r(2) - coeff*r(3); 365 366 coeff = lhs(4,3); 367 lhs(4,4)= lhs(4,4) - coeff*lhs(3,4); 368 lhs(4,5)= lhs(4,5) - coeff*lhs(3,5); 369 c(4,1) = c(4,1) - coeff*c(3,1); 370 c(4,2) = c(4,2) - coeff*c(3,2); 371 c(4,3) = c(4,3) - coeff*c(3,3); 372 c(4,4) = c(4,4) - coeff*c(3,4); 373 c(4,5) = c(4,5) - coeff*c(3,5); 374 r(4) = r(4) - coeff*r(3); 375 376 coeff = lhs(5,3); 377 lhs(5,4)= lhs(5,4) - coeff*lhs(3,4); 378 lhs(5,5)= lhs(5,5) - coeff*lhs(3,5); 379 c(5,1) = c(5,1) - coeff*c(3,1); 380 c(5,2) = c(5,2) - coeff*c(3,2); 381 c(5,3) = c(5,3) - coeff*c(3,3); 382 c(5,4) = c(5,4) - coeff*c(3,4); 383 c(5,5) = c(5,5) - coeff*c(3,5); 384 r(5) = r(5) - coeff*r(3); 385 386 387 pivot = 1.00e0/lhs(4,4); 388 lhs(4,5) = lhs(4,5)*pivot; 389 c(4,1) = c(4,1)*pivot; 390 c(4,2) = c(4,2)*pivot; 391 c(4,3) = c(4,3)*pivot; 392 c(4,4) = c(4,4)*pivot; 393 c(4,5) = c(4,5)*pivot; 394 r(4) = r(4) *pivot; 395 396 coeff = lhs(1,4); 397 lhs(1,5)= lhs(1,5) - coeff*lhs(4,5); 398 c(1,1) = c(1,1) - coeff*c(4,1); 399 c(1,2) = c(1,2) - coeff*c(4,2); 400 c(1,3) = c(1,3) - coeff*c(4,3); 401 c(1,4) = c(1,4) - coeff*c(4,4); 402 c(1,5) = c(1,5) - coeff*c(4,5); 403 r(1) = r(1) - coeff*r(4); 404 405 coeff = lhs(2,4); 406 lhs(2,5)= lhs(2,5) - coeff*lhs(4,5); 407 c(2,1) = c(2,1) - coeff*c(4,1); 408 c(2,2) = c(2,2) - coeff*c(4,2); 409 c(2,3) = c(2,3) - coeff*c(4,3); 410 c(2,4) = c(2,4) - coeff*c(4,4); 411 c(2,5) = c(2,5) - coeff*c(4,5); 412 r(2) = r(2) - coeff*r(4); 413 414 coeff = lhs(3,4); 415 lhs(3,5)= lhs(3,5) - coeff*lhs(4,5); 416 c(3,1) = c(3,1) - coeff*c(4,1); 417 c(3,2) = c(3,2) - coeff*c(4,2); 418 c(3,3) = c(3,3) - coeff*c(4,3); 419 c(3,4) = c(3,4) - coeff*c(4,4); 420 c(3,5) = c(3,5) - coeff*c(4,5); 421 r(3) = r(3) - coeff*r(4); 422 423 coeff = lhs(5,4); 424 lhs(5,5)= lhs(5,5) - coeff*lhs(4,5); 425 c(5,1) = c(5,1) - coeff*c(4,1); 426 c(5,2) = c(5,2) - coeff*c(4,2); 427 c(5,3) = c(5,3) - coeff*c(4,3); 428 c(5,4) = c(5,4) - coeff*c(4,4); 429 c(5,5) = c(5,5) - coeff*c(4,5); 430 r(5) = r(5) - coeff*r(4); 431 432 433 pivot = 1.00e0/lhs(5,5); 434 c(5,1) = c(5,1)*pivot; 435 c(5,2) = c(5,2)*pivot; 436 c(5,3) = c(5,3)*pivot; 437 c(5,4) = c(5,4)*pivot; 438 c(5,5) = c(5,5)*pivot; 439 r(5) = r(5) *pivot; 440 441 coeff = lhs(1,5); 442 c(1,1) = c(1,1) - coeff*c(5,1); 443 c(1,2) = c(1,2) - coeff*c(5,2); 444 c(1,3) = c(1,3) - coeff*c(5,3); 445 c(1,4) = c(1,4) - coeff*c(5,4); 446 c(1,5) = c(1,5) - coeff*c(5,5); 447 r(1) = r(1) - coeff*r(5); 448 449 coeff = lhs(2,5); 450 c(2,1) = c(2,1) - coeff*c(5,1); 451 c(2,2) = c(2,2) - coeff*c(5,2); 452 c(2,3) = c(2,3) - coeff*c(5,3); 453 c(2,4) = c(2,4) - coeff*c(5,4); 454 c(2,5) = c(2,5) - coeff*c(5,5); 455 r(2) = r(2) - coeff*r(5); 456 457 coeff = lhs(3,5); 458 c(3,1) = c(3,1) - coeff*c(5,1); 459 c(3,2) = c(3,2) - coeff*c(5,2); 460 c(3,3) = c(3,3) - coeff*c(5,3); 461 c(3,4) = c(3,4) - coeff*c(5,4); 462 c(3,5) = c(3,5) - coeff*c(5,5); 463 r(3) = r(3) - coeff*r(5); 464 465 coeff = lhs(4,5); 466 c(4,1) = c(4,1) - coeff*c(5,1); 467 c(4,2) = c(4,2) - coeff*c(5,2); 468 c(4,3) = c(4,3) - coeff*c(5,3); 469 c(4,4) = c(4,4) - coeff*c(5,4); 470 c(4,5) = c(4,5) - coeff*c(5,5); 471 r(4) = r(4) - coeff*r(5); 472 473 474 return; 475} 476 477 478 479//--------------------------------------------------------------------- 480//--------------------------------------------------------------------- 481 482void binvrhs( double lhs[],double r[] ) { 483 484//--------------------------------------------------------------------- 485//--------------------------------------------------------------------- 486 487//--------------------------------------------------------------------- 488// 489//--------------------------------------------------------------------- 490 491 double pivot, coeff; 492 493//--------------------------------------------------------------------- 494// 495//--------------------------------------------------------------------- 496 497 498 pivot = 1.00e0/lhs(1,1); 499 lhs(1,2) = lhs(1,2)*pivot; 500 lhs(1,3) = lhs(1,3)*pivot; 501 lhs(1,4) = lhs(1,4)*pivot; 502 lhs(1,5) = lhs(1,5)*pivot; 503 r(1) = r(1) *pivot; 504 505 coeff = lhs(2,1); 506 lhs(2,2)= lhs(2,2) - coeff*lhs(1,2); 507 lhs(2,3)= lhs(2,3) - coeff*lhs(1,3); 508 lhs(2,4)= lhs(2,4) - coeff*lhs(1,4); 509 lhs(2,5)= lhs(2,5) - coeff*lhs(1,5); 510 r(2) = r(2) - coeff*r(1); 511 512 coeff = lhs(3,1); 513 lhs(3,2)= lhs(3,2) - coeff*lhs(1,2); 514 lhs(3,3)= lhs(3,3) - coeff*lhs(1,3); 515 lhs(3,4)= lhs(3,4) - coeff*lhs(1,4); 516 lhs(3,5)= lhs(3,5) - coeff*lhs(1,5); 517 r(3) = r(3) - coeff*r(1); 518 519 coeff = lhs(4,1); 520 lhs(4,2)= lhs(4,2) - coeff*lhs(1,2); 521 lhs(4,3)= lhs(4,3) - coeff*lhs(1,3); 522 lhs(4,4)= lhs(4,4) - coeff*lhs(1,4); 523 lhs(4,5)= lhs(4,5) - coeff*lhs(1,5); 524 r(4) = r(4) - coeff*r(1); 525 526 coeff = lhs(5,1); 527 lhs(5,2)= lhs(5,2) - coeff*lhs(1,2); 528 lhs(5,3)= lhs(5,3) - coeff*lhs(1,3); 529 lhs(5,4)= lhs(5,4) - coeff*lhs(1,4); 530 lhs(5,5)= lhs(5,5) - coeff*lhs(1,5); 531 r(5) = r(5) - coeff*r(1); 532 533 534 pivot = 1.00e0/lhs(2,2); 535 lhs(2,3) = lhs(2,3)*pivot; 536 lhs(2,4) = lhs(2,4)*pivot; 537 lhs(2,5) = lhs(2,5)*pivot; 538 r(2) = r(2) *pivot; 539 540 coeff = lhs(1,2); 541 lhs(1,3)= lhs(1,3) - coeff*lhs(2,3); 542 lhs(1,4)= lhs(1,4) - coeff*lhs(2,4); 543 lhs(1,5)= lhs(1,5) - coeff*lhs(2,5); 544 r(1) = r(1) - coeff*r(2); 545 546 coeff = lhs(3,2); 547 lhs(3,3)= lhs(3,3) - coeff*lhs(2,3); 548 lhs(3,4)= lhs(3,4) - coeff*lhs(2,4); 549 lhs(3,5)= lhs(3,5) - coeff*lhs(2,5); 550 r(3) = r(3) - coeff*r(2); 551 552 coeff = lhs(4,2); 553 lhs(4,3)= lhs(4,3) - coeff*lhs(2,3); 554 lhs(4,4)= lhs(4,4) - coeff*lhs(2,4); 555 lhs(4,5)= lhs(4,5) - coeff*lhs(2,5); 556 r(4) = r(4) - coeff*r(2); 557 558 coeff = lhs(5,2); 559 lhs(5,3)= lhs(5,3) - coeff*lhs(2,3); 560 lhs(5,4)= lhs(5,4) - coeff*lhs(2,4); 561 lhs(5,5)= lhs(5,5) - coeff*lhs(2,5); 562 r(5) = r(5) - coeff*r(2); 563 564 565 pivot = 1.00e0/lhs(3,3); 566 lhs(3,4) = lhs(3,4)*pivot; 567 lhs(3,5) = lhs(3,5)*pivot; 568 r(3) = r(3) *pivot; 569 570 coeff = lhs(1,3); 571 lhs(1,4)= lhs(1,4) - coeff*lhs(3,4); 572 lhs(1,5)= lhs(1,5) - coeff*lhs(3,5); 573 r(1) = r(1) - coeff*r(3); 574 575 coeff = lhs(2,3); 576 lhs(2,4)= lhs(2,4) - coeff*lhs(3,4); 577 lhs(2,5)= lhs(2,5) - coeff*lhs(3,5); 578 r(2) = r(2) - coeff*r(3); 579 580 coeff = lhs(4,3); 581 lhs(4,4)= lhs(4,4) - coeff*lhs(3,4); 582 lhs(4,5)= lhs(4,5) - coeff*lhs(3,5); 583 r(4) = r(4) - coeff*r(3); 584 585 coeff = lhs(5,3); 586 lhs(5,4)= lhs(5,4) - coeff*lhs(3,4); 587 lhs(5,5)= lhs(5,5) - coeff*lhs(3,5); 588 r(5) = r(5) - coeff*r(3); 589 590 591 pivot = 1.00e0/lhs(4,4); 592 lhs(4,5) = lhs(4,5)*pivot; 593 r(4) = r(4) *pivot; 594 595 coeff = lhs(1,4); 596 lhs(1,5)= lhs(1,5) - coeff*lhs(4,5); 597 r(1) = r(1) - coeff*r(4); 598 599 coeff = lhs(2,4); 600 lhs(2,5)= lhs(2,5) - coeff*lhs(4,5); 601 r(2) = r(2) - coeff*r(4); 602 603 coeff = lhs(3,4); 604 lhs(3,5)= lhs(3,5) - coeff*lhs(4,5); 605 r(3) = r(3) - coeff*r(4); 606 607 coeff = lhs(5,4); 608 lhs(5,5)= lhs(5,5) - coeff*lhs(4,5); 609 r(5) = r(5) - coeff*r(4); 610 611 612 pivot = 1.00e0/lhs(5,5); 613 r(5) = r(5) *pivot; 614 615 coeff = lhs(1,5); 616 r(1) = r(1) - coeff*r(5); 617 618 coeff = lhs(2,5); 619 r(2) = r(2) - coeff*r(5); 620 621 coeff = lhs(3,5); 622 r(3) = r(3) - coeff*r(5); 623 624 coeff = lhs(4,5); 625 r(4) = r(4) - coeff*r(5); 626 627 628 return; 629} 630 631 632 633