1# Author : Fabian Jakobs 2r"""BLAS - Basic Linear Algebra Subprograms 3 4GSL provides dense vector and matrix objects, based on the relevant built-in 5types. The library provides an interface to the BLAS operations which apply to 6these objects. 7 8PyGSL only provides the functions working on "native" Python datatypes, i.e. 9double and complex_double. 10 11Functions with the postfix "_cr" are functions that support call by reference. 12This is faster than the original version but may confuse real Python 13programmers. Use this, if speed matters!! 14 15Functions that are naturally done using functions of the underlying numerical 16package are left out here. 17""" 18 19from . import gslwrap 20from . import _gslwrap 21import pygsl._numobj as Numeric 22import copy 23from pygsl import array_typed_copy, get_typecode, Float, Complex 24# 25# constants 26# 27CblasRowMajor = gslwrap.CblasRowMajor 28CblasColMajor = gslwrap.CblasColMajor 29CblasNoTrans = gslwrap.CblasNoTrans 30CblasTrans = gslwrap.CblasTrans 31CblasConjTrans = gslwrap.CblasConjTrans 32CblasUpper = gslwrap.CblasUpper 33CblasLower = gslwrap.CblasLower 34CblasNonUnit = gslwrap.CblasNonUnit 35CblasUnit = gslwrap.CblasUnit 36CblasLeft = gslwrap.CblasLeft 37CblasRight = gslwrap.CblasRight 38gsl_blas_sdsdot = gslwrap.gsl_blas_sdsdot 39 40# 41# Level 1 42# 43 44def ddot(x, y): 45 """This function computes the scalar product \M{x^T y} for the vectors x and y, 46 returning the result. 47 """ 48 return _gslwrap.gsl_blas_ddot(x, y)#[1] 49 50 51def zdotu(x, y): 52 """This function computes the complex scalar product \M{x^T y} for the 53 vectors x and y, returning the result. 54 """ 55 return _gslwrap.gsl_blas_zdotu(x, y, 1j)#[1] 56 57 58def zdotc(x, y): 59 """This function computes the complex conjugate scalar product \M{x^H y} for the 60 vectors x and y, returning the result. 61 """ 62 return _gslwrap.gsl_blas_zdotc(x, y, 1j)#[1] 63 64 65def dnrm2(x): 66 """This function computes the Euclidean norm \M{||x||_2 = \sqrt {\sum x_i^2}} 67 of the vector x. 68 """ 69 return _gslwrap.gsl_blas_dnrm2(x) 70 71 72def dznrm2(x): 73 """This function computes the Euclidean norm of the complex vector x, 74 \M{||x||_2 = \sqrt {\sum (\Re(x_i)^2 + \Im(x_i)^2)}}. 75 """ 76 return _gslwrap.gsl_blas_dznrm2(x) 77 78 79def dasum(x): 80 """This function computes the absolute sum \M{\sum |x_i|} of the elements 81 of the vector x. 82 """ 83 return _gslwrap.gsl_blas_dasum(x) 84 85 86def dzasum(x): 87 """This function computes the absolute sum \M{\sum |\Re(x_i)| + |\Im(x_i)|} 88 of the elements of the vector x. 89 """ 90 return _gslwrap.gsl_blas_dzasum(x) 91 92 93def idamax(x): 94 """This function returns the index of the largest element of the vector x. 95 The largest element is determined by its absolute magnitude. If the 96 largest value occurs several times then the index of the first occurrence 97 is returned. 98 """ 99 return _gslwrap.gsl_blas_idamax(x) 100 101 102def izamax(x): 103 """This function returns the index of the largest element of the vector x. 104 The largest element is determined by the sum of the magnitudes of the 105 real and imaginary parts \M{|\Re(x_i)| + |\Im(x_i)|}. If the largest value 106 occurs several times then the index of the first occurrence is returned. 107 """ 108 return _gslwrap.gsl_blas_izamax(x) 109 110 111def daxpy(alpha, x, y): 112 """This function computes the sum \M{y = S{alpha} x + y} for the vectors x and y. 113 """ 114 yn = array_typed_copy(y, Float) 115 _gslwrap.gsl_blas_daxpy(alpha, x, yn) 116 return yn 117 118 119def daxpy_cr(alpha, x, y_CR): 120 """This function computes the sum \M{y = S{alpha} x + y} for the vectors x and y. 121 """ 122 _gslwrap.gsl_blas_daxpy(alpha, x, y_CR) 123 124 125def zaxpy(alpha, x, y): 126 """This function computes the sum \M{y = S{alpha} x + y} for the vectors x and y. 127 """ 128 yn = array_typed_copy(y, Complex) 129 _gslwrap.gsl_blas_zaxpy(alpha, x, yn) 130 return yn 131 132 133def zaxpy_cr(alpha, x, y_CR): 134 """This function computes the sum \M{y = S{alpha} x + y} for the vectors x and y. 135 """ 136 _gslwrap.gsl_blas_zaxpy(alpha, x, y_CR) 137 138 139def drot(x, y, c, s): 140 """This function applies a Givens rotation \M{(x', y') = (c x + s y, -s x + c y)} 141 to the vectors x, y. 142 """ 143 xn = array_typed_copy(x, Float) 144 yn = array_typed_copy(y, Float) 145 _gslwrap.gsl_blas_drot(xn, yn, c, s) 146 return (xn, yn) 147 148 149def drot_cr(x_CR, y_CR, c, s): 150 """This function applies a Givens rotation \M{(x', y') = (c x + s y, -s x + c y)} 151 to the vectors x, y. 152 """ 153 _gslwrap.gsl_blas_drot(x_CR, y_CR, c, s) 154 155 156# 157# Level 2 158# 159 160def dgemv(alpha, a, x, beta, y, TransA=CblasNoTrans): 161 """This function computes the matrix-vector product and 162 sum \M{y = S{alpha} op(A) x + S{beta} y}, where op(A) = \M{A, A^T, A^H} for 163 TransA = CblasNoTrans, CblasTrans, CblasConjTrans. 164 """ 165 yn = array_typed_copy(y, Float) 166 _gslwrap.gsl_blas_dgemv(TransA, alpha, a, x, beta, yn) 167 return yn 168 169 170def zgemv(alpha, a, x, beta, y, TransA=CblasNoTrans): 171 """This function computes the matrix-vector product and 172 sum \M{y = S{alpha} op(A) x + S{beta} y}, where \M{op(A) = A, A^T, A^H} for 173 TransA = CblasNoTrans, CblasTrans, CblasConjTrans. 174 """ 175 yn = array_typed_copy(y, Complex) 176 _gslwrap.gsl_blas_zgemv(TransA, alpha, a, x, beta, yn) 177 return yn 178 179 180def dtrmv(A, 181 x, 182 Uplo=CblasLower, 183 TransA=CblasNoTrans, 184 Diag=CblasNonUnit): 185 """This function computes the matrix-vector product 186 returns x' 187 188 This function computes the matrix-vector product and 189 x' = op(A) x for the triangular 190 matrix A, where op(A) = A, A^T, A^H for TransA = CblasNoTrans, 191 CblasTrans, CblasConjTrans. When Uplo is CblasUpper then the upper 192 triangle of A is used, and when Uplo is CblasLower then the lower 193 triangle of A is used. If Diag is CblasNonUnit then the diagonal of 194 the matrix is used, but if Diag is CblasUnit then the diagonal 195 elements of the matrix A are taken as unity and are not referenced. 196 """ 197 xn = array_typed_copy(x, Float) 198 199 _gslwrap.gsl_blas_dtrmv(Uplo, TransA, Diag, A, xn) 200 return xn 201 202 203def ztrmv(A, 204 x, 205 Uplo=CblasLower, 206 TransA=CblasNoTrans, 207 Diag=CblasNonUnit): 208 """ 209 returns x' 210 211 This function computes the matrix-vector product and 212 x' = op(A) x for the triangular 213 matrix A, where op(A) = A, A^T, A^H for TransA = CblasNoTrans, 214 CblasTrans, CblasConjTrans. When Uplo is CblasUpper then the upper 215 triangle of A is used, and when Uplo is CblasLower then the lower 216 triangle of A is used. If Diag is CblasNonUnit then the diagonal of 217 the matrix is used, but if Diag is CblasUnit then the diagonal 218 elements of the matrix A are taken as unity and are not referenced. 219 """ 220 xn = array_typed_copy(x) 221 _gslwrap.gsl_blas_ztrmv(Uplo, TransA, Diag, A, xn) 222 return xn 223 224 225def dtrsv(A, 226 x, 227 Uplo=CblasLower, 228 TransA=CblasNoTrans, 229 Diag=CblasNonUnit): 230 """ 231 232 This function computes inv(op(A)) x for x, where op(A) = A, A^T, A^H 233 for TransA = CblasNoTrans, CblasTrans, CblasConjTrans. When Uplo is 234 CblasUpper then the upper triangle of A is used, and when Uplo is 235 CblasLower then the lower triangle of A is used. If Diag is 236 CblasNonUnit then the diagonal of the matrix is used, but if Diag 237 is CblasUnit then the diagonal elements of the matrix A are taken 238 as unity and are not referenced. 239 240 241 242 Returns: 243 x: 244 """ 245 xn = array_typed_copy(x) 246 _gslwrap.gsl_blas_dtrsv(Uplo, TransA, Diag, A, xn) 247 return xn 248 249 250def ztrsv(A, 251 x, 252 Uplo=CblasLower, 253 TransA=CblasNoTrans, 254 Diag=CblasNonUnit): 255 """ 256 returns x' 257 258 This function computes inv(op(A)) x for x, where op(A) = A, A^T, A^H 259 for TransA = CblasNoTrans, CblasTrans, CblasConjTrans. When Uplo is 260 CblasUpper then the upper triangle of A is used, and when Uplo is 261 CblasLower then the lower triangle of A is used. If Diag is 262 CblasNonUnit then the diagonal of the matrix is used, but if Diag 263 is CblasUnit then the diagonal elements of the matrix A are taken 264 as unity and are not referenced. 265 """ 266 xn = array_typed_copy(x) 267 _gslwrap.gsl_blas_ztrsv(Uplo, TransA, Diag, A, xn) 268 return xn 269 270 271def dsymv(alpha, A, X, beta, Y, Uplo=CblasLower): 272 """ 273 returns y' 274 275 This function computes the matrix-vector product and 276 sum \M{y' = S{alpha} A x + S{beta} y} for the symmetric matrix A. Since the 277 matrix A is symmetric only its upper half or lower half need to be 278 stored. When Uplo is CblasUpper then the upper triangle and diagonal 279 of A are used, and when Uplo is CblasLower then the lower triangle 280 and diagonal of A are used. 281 """ 282 yn = array_typed_copy(Y, Float) 283 _gslwrap.gsl_blas_dsymv(Uplo, alpha, A, X, beta, yn) 284 return yn 285 286 287def zhemv(alpha, A, X, beta, Y, Uplo=CblasLower): 288 """ 289 returns y' 290 291 This function computes the matrix-vector product and 292 sum \M{y' = S{alpha} A x + S{beta} y} for the hermitian matrix A. Since the 293 matrix A is hermitian only its upper half or lower half need to be 294 stored. When Uplo is CblasUpper then the upper triangle and diagonal 295 of A are used, and when Uplo is CblasLower then the lower triangle 296 and diagonal of A are used. The imaginary elements of the diagonal 297 are automatically assumed to be zero and are not referenced. 298 """ 299 yn = array_typed_copy(Y, get_typecode(A)) 300 _gslwrap.gsl_blas_zhemv(Uplo, alpha, A, X, beta, yn) 301 return yn 302 303 304def dger(alpha, X, Y, A): 305 """ 306 returns A' 307 308 This function computes the rank-1 update \M{A' = S{alpha} x y^T + A} of 309 the matrix A. 310 """ 311 an = array_typed_copy(A) 312 _gslwrap.gsl_blas_dger(alpha, X, Y, an) 313 return an 314 315 316def zgeru(alpha, X, Y, A): 317 """ 318 returns A' 319 320 This function computes the rank-1 update \M{A' = S{alpha} x y^T + A} of 321 the matrix A. 322 """ 323 an = array_typed_copy(A) 324 _gslwrap.gsl_blas_zgeru(alpha, X, Y, an) 325 return an 326 327 328def zgerc(alpha, X, Y, A): 329 """ 330 returns A' 331 332 This function computes the conjugate rank-1 update 333 \M{A = S{alpha} x y^H + A} of the matrix A. 334 """ 335 an = array_typed_copy(A) 336 _gslwrap.gsl_blas_zgerc(alpha, X, Y, an) 337 return an 338 339 340def dsyr(alpha, X, A, Uplo=CblasLower): 341 """ 342 returns A' 343 344 This function computes the symmetric rank-1 update 345 \M{A' = S{alpha} x x^T + A} of the symmetric matrix A. Since the matrix A 346 is symmetric only its upper half or lower half need to be stored. 347 When Uplo is CblasUpper then the upper triangle and diagonal of A 348 are used, and when Uplo is CblasLower then the lower triangle and 349 diagonal of A are used. 350 """ 351 an = array_typed_copy(A) 352 _gslwrap.gsl_blas_dsyr(Uplo, alpha, X, an) 353 return an 354 355 356def zher(alpha, X, A, Uplo=CblasLower): 357 """ 358 returns A' 359 360 This function computes the hermitian rank-1 update 361 \M{A' = S{alpha} x x^H + A} of the hermitian matrix A. Since the matrix A 362 is hermitian only its upper half or lower half need to be stored. 363 When Uplo is CblasUpper then the upper triangle and diagonal of A are 364 used, and when Uplo is CblasLower then the lower triangle and diagonal 365 of A are used. The imaginary elements of the diagonal are automatically 366 set to zero. 367 """ 368 an = array_typed_copy(A) 369 _gslwrap.gsl_blas_zher(Uplo, alpha, X, an) 370 return an 371 372 373def dsyr2(alpha, X, Y, A, Uplo=CblasLower): 374 """ 375 returns A' 376 377 This function computes the symmetric rank-2 update 378 \M{A' = S{alpha} x y^T + S{alpha} y x^T + A} of the symmetric matrix A. 379 Since the matrix A is symmetric only its upper half or lower half 380 need to be stored. When Uplo is CblasUpper then the upper triangle 381 and diagonal of A are used, and when Uplo is CblasLower then the 382 lower triangle and diagonal of A are used. 383 """ 384 an = array_typed_copy(A) 385 _gslwrap.gsl_blas_dsyr2(Uplo, alpha, X, Y, an) 386 return an 387 388 389def zher2(alpha, X, Y, A, Uplo=CblasLower): 390 """ 391 returns A' 392 393 This function computes the hermitian rank-2 update 394 \M{A' = S{alpha} x y^H + S{alpha}^* y x^H A} of the hermitian matrix A. 395 Since the matrix A is hermitian only its upper half or lower half 396 need to be stored. When Uplo is CblasUpper then the upper triangle 397 and diagonal of A are used, and when Uplo is CblasLower then the 398 lower triangle and diagonal of A are used. The imaginary elements 399 of the diagonal are automatically set to zero. 400 """ 401 an = array_typed_copy(A) 402 _gslwrap.gsl_blas_zher2(Uplo, alpha, X, Y, an) 403 return an 404 405 406# 407# Level 3 408# 409 410def dgemm(alpha, 411 A, B, 412 beta, C, 413 TransA=CblasNoTrans, 414 TransB=CblasNoTrans): 415 """ 416 returns C' 417 418 This function computes the matrix-matrix product and sum 419 \M{C' = S{alpha} op(A) op(B) + S{beta} C} where op(A) = A, A^T, A^H for 420 TransA = CblasNoTrans, CblasTrans, CblasConjTrans and similarly 421 for the parameter TransB. 422 """ 423 cn = array_typed_copy(C) 424 _gslwrap.gsl_blas_dgemm(TransA, TransB, alpha, A, B, beta, cn) 425 return cn 426 427 428def zgemm(alpha, 429 A, B, 430 beta, C, 431 TransA=CblasNoTrans, 432 TransB=CblasNoTrans): 433 """ 434 returns C' 435 436 This function computes the matrix-matrix product and sum 437 \M{C' = S{alpha} op(A) op(B) + S{beta} C} where op(A) = A, A^T, A^H for 438 TransA = CblasNoTrans, CblasTrans, CblasConjTrans and similarly 439 for the parameter TransB. 440 """ 441 cn = array_typed_copy(C) 442 _gslwrap.gsl_blas_zgemm(TransA, TransB, alpha, A, B, beta, cn) 443 return cn 444 445 446def dsymm(alpha, A, B, beta, C, 447 Side=CblasLeft, 448 Uplo=CblasLower): 449 """ 450 returns C' 451 452 This function computes the matrix-matrix product and 453 sum \M{C' = S{alpha} A B + S{beta} C} for Side is CblasLeft and 454 \M{C' = S{alpha} B A + S{beta} C} for Side is CblasRight, where the matrix A 455 is symmetric. When Uplo is CblasUpper then the upper triangle and 456 diagonal of A are used, and when Uplo is CblasLower then the lower 457 triangle and diagonal of A are used. 458 """ 459 cn = array_typed_copy(C) 460 _gslwrap.gsl_blas_dsymm(Side, Uplo, alpha, A, B, beta, cn) 461 return cn 462 463 464def zsymm(alpha, A, B, beta, C, 465 Side=CblasLeft, 466 Uplo=CblasLower): 467 """ 468 returns C' 469 470 This function computes the matrix-matrix product and 471 sum \M{C' = S{alpha} A B + S{beta} C} for Side is CblasLeft and 472 \M{C' = S{alpha} B A + S{beta} C} for Side is CblasRight, where the matrix A 473 is symmetric. When Uplo is CblasUpper then the upper triangle and 474 diagonal of A are used, and when Uplo is CblasLower then the lower 475 triangle and diagonal of A are used. 476 """ 477 cn = array_typed_copy(C) 478 _gslwrap.gsl_blas_zsymm(Side, Uplo, alpha, A, B, beta, cn) 479 return cn 480 481 482def zhemm(alpha, A, B, beta, C, 483 Side=CblasLeft, 484 Uplo=CblasLower): 485 """ 486 returns C' 487 488 This function computes the matrix-matrix product and 489 sum \M{C' = S{alpha} A B + S{beta} C} for Side is CblasLeft and 490 \M{C' = S{alpha} B A + S{beta} C} for Side is CblasRight, where the matrix A 491 is hermitian. When Uplo is CblasUpper then the upper triangle and 492 diagonal of A are used, and when Uplo is CblasLower then the lower 493 triangle and diagonal of A are used. The imaginary elements of the 494 diagonal are automatically set to zero. 495 """ 496 cn = array_typed_copy(C) 497 _gslwrap.gsl_blas_zhemm(Side, Uplo, alpha, A, B, beta, cn) 498 return cn 499 500 501def dtrmm(alpha, A, B, 502 Side=CblasLeft, 503 Uplo=CblasLower, 504 TransA=CblasNoTrans, 505 Diag=CblasNonUnit): 506 """ 507 returns B' 508 509 This function computes the matrix-matrix product 510 \M{B' = S{alpha} op(A) B} for Side is CblasLeft and 511 \M{B' = S{alpha} B op(A)} for Side is CblasRight. The matrix A is 512 triangular and op(A) = A, A^T, A^H for TransA = CblasNoTrans, 513 CblasTrans, CblasConjTrans When Uplo is CblasUpper then the upper 514 triangle of A is used, and when Uplo is CblasLower then the lower 515 triangle of A is used. If Diag is CblasNonUnit then the diagonal 516 of A is used, but if Diag is CblasUnit then the diagonal elements 517 of the matrix A are taken as unity and are not referenced. 518 """ 519 bn = array_typed_copy(B) 520 _gslwrap.gsl_blas_dtrmm(Side, Uplo, TransA, Diag, alpha, A, bn) 521 return bn 522 523 524def ztrmm(alpha, A, B, 525 Side=CblasLeft, 526 Uplo=CblasLower, 527 TransA=CblasNoTrans, 528 Diag=CblasNonUnit): 529 """ 530 returns B' 531 532 This function computes the matrix-matrix product 533 \M{B' = S{alpha} op(A) B} for Side is CblasLeft and 534 \M{B' = S{alpha} B op(A)} for Side is CblasRight. The matrix A is 535 triangular and op(A) = A, A^T, A^H for TransA = CblasNoTrans, 536 CblasTrans, CblasConjTrans When Uplo is CblasUpper then the upper 537 triangle of A is used, and when Uplo is CblasLower then the lower 538 triangle of A is used. If Diag is CblasNonUnit then the diagonal 539 of A is used, but if Diag is CblasUnit then the diagonal elements 540 of the matrix A are taken as unity and are not referenced. 541 """ 542 bn = array_typed_copy(B) 543 _gslwrap.gsl_blas_ztrmm(Side, Uplo, TransA, Diag, alpha, A, bn) 544 return bn 545 546 547def dtrsm(alpha, A, B, 548 Side=CblasLeft, 549 Uplo=CblasLower, 550 TransA=CblasNoTrans, 551 Diag=CblasNonUnit): 552 """ 553 returns B' 554 555 This function computes the matrix-matrix product 556 \M{B' = S{alpha} op(inv(A)) B} for Side is CblasLeft and 557 \M{B' = S{alpha} B op(inv(A))} for Side is CblasRight. The matrix A is 558 triangular and op(A) = A, A^T, A^H for TransA = CblasNoTrans, 559 CblasTrans, CblasConjTrans When Uplo is CblasUpper then the upper 560 triangle of A is used, and when Uplo is CblasLower then the lower 561 triangle of A is used. If Diag is CblasNonUnit then the diagonal 562 of A is used, but if Diag is CblasUnit then the diagonal elements 563 of the matrix A are taken as unity and are not referenced. 564 """ 565 bn = array_typed_copy(B) 566 _gslwrap.gsl_blas_dtrsm(Side, Uplo, TransA, Diag, alpha, A, bn) 567 return bn 568 569 570def ztrsm(alpha, A, B, 571 Side=CblasLeft, 572 Uplo=CblasLower, 573 TransA=CblasNoTrans, 574 Diag=CblasNonUnit): 575 """ 576 Returns: 577 :math:`B'` 578 579 This function computes the matrix-matrix product 580 \M{B' = S{alpha} op(inv(A)) B} for Side is CblasLeft and 581 \M{ B' = S{alpha} B op(inv(A))} for Side is CblasRight. The matrix A is 582 triangular and op(A) = A, A^T, A^H for TransA = CblasNoTrans, 583 CblasTrans, CblasConjTrans When Uplo is CblasUpper then the upper 584 triangle of A is used, and when Uplo is CblasLower then the lower 585 triangle of A is used. If Diag is CblasNonUnit then the diagonal 586 of A is used, but if Diag is CblasUnit then the diagonal elements 587 of the matrix A are taken as unity and are not referenced. 588 """ 589 bn = array_typed_copy(B) 590 _gslwrap.gsl_blas_ztrsm(Side, Uplo, TransA, Diag, alpha, A, bn) 591 return bn 592 593 594def dsyrk(alpha, A, beta, C, 595 Uplo=CblasLower, 596 Trans=CblasNoTrans): 597 """ 598 returns C' 599 600 This function computes a rank-k update of the symmetric matrix C, 601 \M{C' = S{alpha} A A^T + S{beta} C} when Trans is CblasNoTrans and 602 \M{C' = S{alpha} A^T A + S{beta} C} when Trans is CblasTrans. Since the matrix 603 C is symmetric only its upper half or lower half need to be stored. 604 When Uplo is CblasUpper then the upper triangle and diagonal of C are 605 used, and when Uplo is CblasLower then the lower triangle and diagonal 606 of C are used. 607 """ 608 cn = array_typed_copy(C) 609 _gslwrap.gsl_blas_dsyrk(Uplo, Trans, alpha, A, beta, cn) 610 return cn 611 612 613def zsyrk(alpha, A, beta, C, 614 Uplo=CblasLower, 615 Trans=CblasNoTrans): 616 """ 617 returns C' 618 619 This function computes a rank-k update of the symmetric matrix C, 620 \M{C' = S{alpha} A A^T + S{beta} C} when Trans is CblasNoTrans and 621 \M{C' = S{alpha} A^T A + S{beta} C} when Trans is CblasTrans. Since the matrix 622 C is symmetric only its upper half or lower half need to be stored. 623 When Uplo is CblasUpper then the upper triangle and diagonal of C are 624 used, and when Uplo is CblasLower then the lower triangle and diagonal 625 of C are used. 626 """ 627 cn = array_typed_copy(C) 628 _gslwrap.gsl_blas_zsyrk(Uplo, Trans, alpha, A, beta, cn) 629 return cn 630 631 632def triang2symm(A, 633 Uplo=CblasLower, 634 Diag=CblasNonUnit): 635 """ 636 returns A' 637 638 This function creates an symmetric matrix from a triangular matrix. 639 When Uplo is CblasUpper then the upper triangle and diagonal of A 640 are used, and when Uplo is CblasLower then the lower triangle and 641 diagonal of A are used. If Diag is CblasNonUnit then the diagonal 642 of A is used, but if Diag is CblasUnit then the diagonal elements 643 of the matrix A are taken as unity and are not referenced. 644 """ 645 an = array_typed_copy(A) 646 if Uplo == CblasLower: 647 for i in range(A.shape[0]): 648 an[:i+1,i] = A[i,:i+1] 649 else: 650 for i in range(A.shape[0]): 651 an[i:,i] = A[i,i:] 652 if Diag == CblasUnit: 653 for i in range(an.shape[0]): 654 an[i,i] = 1 655 return an 656 657 658def triang2herm(A, 659 Uplo=CblasLower, 660 Diag=CblasNonUnit): 661 """ 662 returns A' 663 664 This function creates an hermitian matrix from a triangular matrix. 665 When Uplo is CblasUpper then the upper triangle and diagonal of A 666 are used, and when Uplo is CblasLower then the lower triangle and 667 diagonal of A are used. If Diag is CblasNonUnit then the diagonal 668 of A is used, but if Diag is CblasUnit then the diagonal elements 669 of the matrix A are taken as unity and are not referenced. 670 """ 671 an = array_typed_copy(A) 672 if Uplo == CblasLower: 673 for i in range(A.shape[0]): 674 an[:i+1,i] = Numeric.conjugate(A[i,:i+1]) 675 else: 676 for i in range(A.shape[0]): 677 an[i:,i] = Numeric.conjugate(A[i,i:]) 678 if Diag == CblasUnit: 679 for i in range(an.shape[0]): 680 an[i,i] = 1 681 return an 682 683 684def zherk(alpha, A, beta, C, 685 Uplo=CblasLower, 686 Trans=CblasNoTrans): 687 """ 688 returns C' 689 690 This function computes a rank-k update of the hermitian matrix C, 691 \M{C' = S{alpha} A A^H + S{beta} C} when Trans is CblasNoTrans and 692 \M{C' = S{alpha} A^H A + S{beta} C} when Trans is CblasTrans. Since 693 the matrix C is hermitian only its upper half or lower half need to 694 be stored. When Uplo is CblasUpper then the upper triangle and diagonal 695 of C are used, and when Uplo is CblasLower then the lower triangle and 696 diagonal of C are used. The imaginary elements of the diagonal are 697 automatically set to zero. 698 """ 699 cn = array_typed_copy(C) 700 _gslwrap.gsl_blas_zherk(Uplo, Trans, alpha, A, beta, cn) 701 return cn 702 703 704def dsyr2k(alpha, A, B, beta, C, 705 Uplo=CblasLower, 706 Trans=CblasNoTrans): 707 """ 708 returns C' 709 710 This function computes a rank-2k update of the symmetric 711 matrix C, \M{C' = S{alpha} A B^T + S{alpha} B A^T + S{beta} C} when Trans 712 is CblasNoTrans and \M{C' = S{alpha} A^T B + S{alpha} B^T A + S{beta} C} when 713 Trans is CblasTrans. Since the matrix C is symmetric only its upper 714 half or lower half need to be stored. When Uplo is CblasUpper then 715 the upper triangle and diagonal of C are used, and when Uplo is 716 CblasLower then the lower triangle and diagonal of C are used. 717 """ 718 cn = array_typed_copy(C) 719 _gslwrap.gsl_blas_dsyr2k(Uplo, Trans, alpha, A, B, beta, cn) 720 return cn 721 722 723def zsyr2k(alpha, A, B, beta, C, 724 Uplo=CblasLower, 725 Trans=CblasNoTrans): 726 """ 727 returns C' 728 729 This function computes a rank-2k update of the symmetric 730 matrix C, \M{C' = S{alpha} A B^T + S{alpha} B A^T + S{beta} C} when Trans 731 is CblasNoTrans and \M{C' = S{alpha} A^T B + S{alpha} B^T A + S{beta} C} when 732 Trans is CblasTrans. Since the matrix C is symmetric only its upper 733 half or lower half need to be stored. When Uplo is CblasUpper then 734 the upper triangle and diagonal of C are used, and when Uplo is 735 CblasLower then the lower triangle and diagonal of C are used. 736 """ 737 cn = array_typed_copy(C) 738 _gslwrap.gsl_blas_zsyr2k(Uplo, Trans, alpha, A, B, beta, cn) 739 return cn 740 741 742def zher2k(alpha, A, B, beta, C, 743 Uplo=CblasLower, 744 Trans=CblasNoTrans): 745 """ 746 returns C' 747 748 This function computes a rank-2k update of the hermitian matrix C, 749 \M{C' = S{alpha} A B^H + S{alpha}^* B A^H + S{beta} C} when Trans is 750 CblasNoTrans and \M{C' = S{alpha} A^H B + S{alpha}^* B^H A + S{beta} C} when 751 Trans is CblasTrans. Since the matrix C is hermitian only its upper 752 half or lower half need to be stored. When Uplo is CblasUpper then 753 the upper triangle and diagonal of C are used, and when Uplo is 754 CblasLower then the lower triangle and diagonal of C are used. 755 The imaginary elements of the diagonal are automatically set to zero. 756 757 Calls function: 758 """ 759 cn = array_typed_copy(C) 760 _gslwrap.gsl_blas_zher2k(Uplo, Trans, alpha, A, B, beta, cn) 761 return cn 762