1!--------------------------------------------------------------------------------------------------! 2! DFTB+: general package for performing fast atomistic simulations ! 3! Copyright (C) 2006 - 2019 DFTB+ developers group ! 4! ! 5! See the LICENSE file for terms of usage and distribution. ! 6!--------------------------------------------------------------------------------------------------! 7 8#:include 'common.fypp' 9 10!> Contains F90 wrapper functions for some commonly used blas calls needed in the code. The 11!> interface of all BLAS calls must be defined in the module blas. 12module dftbp_blasroutines 13 use dftbp_assert 14 use dftbp_accuracy 15 use dftbp_blas 16 implicit none 17 18 19 !> Rank 1 update of a matrix A := alpha*x*x' + A 20 !> Wrapper for the level 2 blas routine xsyr to perform the rank 1 update of the chosen triangle 21 interface her 22 module procedure her_real 23 module procedure her_cmplx 24 module procedure her_dble 25 module procedure her_dblecmplx 26 end interface her 27 28 29 !> Rank 1 update of a matrix A := alpha*x*y' + A 30 !> Wrapper for the level 2 blas routine xger to perform the rank 1 update of a general matrix 31 interface ger 32 module procedure ger_real 33 module procedure ger_cmplx 34 module procedure ger_dble 35 module procedure ger_dblecmplx 36 end interface ger 37 38 39 !> Symmetric matrix vector multiply y := alpha*A*x + beta*y 40 !> Wrapper for the level 2 blas routine 41 interface hemv 42 module procedure symv_real 43 module procedure symv_dble 44 module procedure hemv_cmplx 45 module procedure hemv_dblecmplx 46 end interface hemv 47 48 49 !> General matrix vector multiply y := alpha*a*x + beta*y 50 !> Wrapper for the level 2 blas routine 51 interface gemv 52 module procedure gemv_real 53 module procedure gemv_dble 54 end interface gemv 55 56 57 !> Banded symmetric matrix vector multiply y := alpha*A*x + beta*y 58 !> Wrapper for the level 2 blas routine 59 interface sbmv 60 module procedure sbmv_real 61 module procedure sbmv_dble 62 end interface sbmv 63 64 65 !> Interface to SYMM routines 66 !> Wrapper for the level 3 blas routine 67 interface symm 68 module procedure symm_real 69 module procedure symm_dble 70 end interface symm 71 72 73 !> Interface to GEMM routines evaluates C := alpha*op( A )*op( B ) + beta*C, where op( X ) is one 74 !> of op( X ) = X or op( X ) = X' 75 76 !> Wrapper for the level 3 blas routine 77 interface gemm 78 module procedure gemm_real 79 module procedure gemm_dble 80 module procedure gemm_cmplx 81 module procedure gemm_dblecmplx 82 end interface gemm 83 84 85 !> Rank k update of a matrix C := alpha*A*A' + beta C 86 !> Wrapper for the level 3 blas routine syrk to perform the rank k update of the chosen triangle 87 !> ofC 88 interface herk 89 module procedure herk_real 90 module procedure herk_cmplx 91 module procedure herk_dble 92 module procedure herk_dblecmplx 93 end interface herk 94 95 96 !> Interface to HEMM routines 97 !> Wrapper for the level 3 blas routine 98 interface hemm 99 module procedure hemm_cmplx 100 module procedure hemm_dblecmplx 101 end interface hemm 102 103contains 104 105 106 !> Real rank 1 update of a symmetric matrix 107 subroutine her_real(a,alpha,x,uplo) 108 109 !> contains the matrix for the update 110 real(rsp), intent(inout) :: a(:,:) 111 112 !> scaling value for the update contribution 113 real(rsp), intent(in) :: alpha 114 115 !> vector of values for the update 116 real(rsp), intent(in) :: x(:) 117 118 !> optional upper, 'U', or lower 'L' triangle, defaults to lower 119 character, intent(in), optional :: uplo 120 121 integer :: n 122 character :: iUplo 123 @:ASSERT(size(a,dim=1) == size(a,dim=2)) 124 @:ASSERT(size(a,dim=1) == size(x)) 125 126 if (present(uplo)) then 127 iuplo = uplo 128 else 129 iuplo = 'l' 130 end if 131 @:ASSERT(iuplo == 'u' .or. iuplo == 'U' .or. iuplo == 'l' .or. iuplo == 'L') 132 n = size(x) 133 call ssyr(iuplo,n,alpha,x,1,a,n) 134 end subroutine her_real 135 136 137 !> Complex rank 1 update of a symmetric matrix 138 subroutine her_cmplx(a,alpha,x,uplo) 139 140 !> contains the matrix for the update 141 complex(rsp), intent(inout) :: a(:,:) 142 143 !> scaling value for the update contribution 144 real(rsp), intent(in) :: alpha 145 146 !> vector of values for the update 147 complex(rsp), intent(in) :: x(:) 148 149 !> optional upper, 'U', or lower 'L' triangle, defaults to lower 150 character, intent(in),optional :: uplo 151 152 integer :: n 153 character :: iuplo 154 @:ASSERT(size(a,dim=1) == size(a,dim=2)) 155 @:ASSERT(size(a,dim=1) == size(x)) 156 157 if (present(uplo)) then 158 iuplo = uplo 159 else 160 iuplo = 'l' 161 end if 162 @:ASSERT(iuplo == 'u' .or. iuplo == 'U' .or. iuplo == 'l' .or. iuplo == 'L') 163 n = size(x) 164 call cher(iuplo,n,alpha,x,1,a,n) 165 end subroutine her_cmplx 166 167 168 !> Double precision rank 1 update of a symmetric matrix 169 subroutine her_dble(a,alpha,x,uplo) 170 171 !> contains the matrix for the update 172 real(rdp), intent(inout) :: a(:,:) 173 174 !> scaling value for the update contribution 175 real(rdp), intent(in) :: alpha 176 177 !> vector of values for the update 178 real(rdp), intent(in) :: x(:) 179 180 !> optional upper, 'U', or lower 'L' triangle, defaults to lower 181 character, intent(in), optional :: uplo 182 183 character :: iuplo 184 integer :: n 185 @:ASSERT(size(a,dim=1) == size(a,dim=2)) 186 @:ASSERT(size(a,dim=1) == size(x)) 187 if (present(uplo)) then 188 iuplo = uplo 189 else 190 iuplo = 'l' 191 end if 192 @:ASSERT(iuplo == 'u' .or. iuplo == 'U' .or. iuplo == 'l' .or. iuplo == 'L') 193 n = size(x) 194 call dsyr(iuplo,n,alpha,x,1,a,n) 195 end subroutine her_dble 196 197 198 !> Double complex rank 1 update of a symmetric matrix 199 subroutine her_dblecmplx(a,alpha,x,uplo) 200 201 !> contains the matrix for the update 202 complex(rdp), intent(inout) :: a(:,:) 203 204 !> scaling value for the update contribution 205 real(rdp), intent(in) :: alpha 206 207 !> vector of values for the update 208 complex(rdp), intent(in) :: x(:) 209 210 !> optional upper, 'U', or lower 'L' triangle, defaults to lower 211 character, intent(in),optional :: uplo 212 213 integer :: n 214 character :: iuplo 215 @:ASSERT(size(a,dim=1) == size(a,dim=2)) 216 @:ASSERT(size(a,dim=1) == size(x)) 217 if (present(uplo)) then 218 iuplo = uplo 219 else 220 iuplo = 'l' 221 end if 222 @:ASSERT(iuplo == 'u' .or. iuplo == 'U' .or. iuplo == 'l' .or. iuplo == 'L') 223 n = size(x) 224 call zher(iuplo,n,alpha,x,1,a,n) 225 end subroutine her_dblecmplx 226 227 228 !> Real rank 1 update of a general matrix 229 subroutine ger_real(a,alpha,x,y) 230 231 !> contains the matrix for the update 232 real(rsp), intent(inout) :: a(:,:) 233 234 !> scaling value for the update contribution 235 real(rsp), intent(in) :: alpha 236 237 !> vector of values for the update 238 real(rsp), intent(in) :: x(:) 239 240 !> vector of values for the update 241 real(rsp), intent(in) :: y(:) 242 243 integer :: n, m 244 @:ASSERT(size(a,dim=1) == size(x)) 245 @:ASSERT(size(a,dim=2) == size(y)) 246 m = size(x) 247 n = size(y) 248 call sger(m,n,alpha,x,1,y,1,a,m) 249 end subroutine ger_real 250 251 252 !> complex rank 1 update of a general matrix 253 subroutine ger_cmplx(a,alpha,x,y) 254 255 !> contains the matrix for the update 256 complex(rsp), intent(inout) :: a(:,:) 257 258 !> scaling value for the update contribution 259 complex(rsp), intent(in) :: alpha 260 261 !> vector of values for the update 262 complex(rsp), intent(in) :: x(:) 263 264 !> vector of values for the update 265 complex(rsp), intent(in) :: y(:) 266 267 integer :: n, m 268 @:ASSERT(size(a,dim=1) == size(x)) 269 @:ASSERT(size(a,dim=2) == size(y)) 270 m = size(x) 271 n = size(y) 272 call cgerc(m,n,alpha,x,1,y,1,a,m) 273 end subroutine ger_cmplx 274 275 276 !> Double precision rank 1 update of a general matrix 277 subroutine ger_dble(a,alpha,x,y) 278 279 !> contains the matrix for the update 280 real(rdp), intent(inout) :: a(:,:) 281 282 !> scaling value for the update contribution 283 real(rdp), intent(in) :: alpha 284 285 !> vector of values for the update 286 real(rdp), intent(in) :: x(:) 287 288 !> vector of values for the update 289 real(rdp), intent(in) :: y(:) 290 291 integer :: n, m 292 @:ASSERT(size(a,dim=1) == size(x)) 293 @:ASSERT(size(a,dim=2) == size(y)) 294 m = size(x) 295 n = size(y) 296 call dger(m,n,alpha,x,1,y,1,a,m) 297 end subroutine ger_dble 298 299 300 !> complex rank 1 update of a general matrix 301 subroutine ger_dblecmplx(a,alpha,x,y) 302 303 !> contains the matrix for the update 304 complex(rdp), intent(inout) :: a(:,:) 305 306 !> scaling value for the update contribution 307 complex(rdp), intent(in) :: alpha 308 309 !> vector of values for the update 310 complex(rdp), intent(in) :: x(:) 311 312 !> vector of values for the update 313 complex(rdp), intent(in) :: y(:) 314 315 integer :: n, m 316 @:ASSERT(size(a,dim=1) == size(x)) 317 @:ASSERT(size(a,dim=2) == size(y)) 318 m = size(x) 319 n = size(y) 320 call zgerc(m,n,alpha,x,1,y,1,a,m) 321 end subroutine ger_dblecmplx 322 323 324 !> real symmetric matrix*vector product 325 subroutine symv_real(y,a,x,uplo,alpha,beta) 326 327 !> vector 328 real(rsp), intent(inout) :: y(:) 329 330 !> symmetric matrix 331 real(rsp), intent(in) :: a(:,:) 332 333 !> vector 334 real(rsp), intent(in) :: x(:) 335 336 !> optional upper, 'U', or lower 'L' triangle (defaults to lower) 337 character, intent(in), optional :: uplo 338 339 !> optional scaling factor (defaults to 1) 340 real(rsp), intent(in), optional :: alpha 341 342 !> optional scaling factor (defaults to 0) 343 real(rsp), intent(in), optional :: beta 344 345 integer :: n 346 character :: iUplo 347 real(rsp) :: iAlpha, iBeta 348 349 if (present(uplo)) then 350 iUplo = uplo 351 else 352 iUplo = 'L' 353 end if 354 if (present(alpha)) then 355 iAlpha = alpha 356 else 357 iAlpha = 1.0_rsp 358 end if 359 if (present(beta)) then 360 iBeta = beta 361 else 362 iBeta = 0.0_rsp 363 end if 364 365 @:ASSERT(size(y) == size(x)) 366 @:ASSERT(size(a,dim=1) == size(a,dim=2)) 367 @:ASSERT(size(a,dim=1) == size(x)) 368 @:ASSERT(iUplo == 'u' .or. iUplo == 'U' .or. iUplo == 'l' .or. iUplo == 'L') 369 n = size(y) 370 call ssymv( iUplo, n, iAlpha, a, n, x, 1, iBeta, y, 1 ) 371 372 end subroutine symv_real 373 374 375 !> real symmetric matrix*vector product 376 subroutine symv_dble(y,a,x,uplo,alpha,beta) 377 378 !> vector 379 real(rdp), intent(inout) :: y(:) 380 381 !> symmetric matrix 382 real(rdp), intent(in) :: a(:,:) 383 384 !> vector 385 real(rdp), intent(in) :: x(:) 386 387 !> optional upper, 'U', or lower 'L' triangle (defaults to lower) 388 character, intent(in), optional :: uplo 389 390 !> optional scaling factor (defaults to 1) 391 real(rdp), intent(in), optional :: alpha 392 393 !> optional scaling factor (defaults to 0) 394 real(rdp), intent(in), optional :: beta 395 396 integer :: n 397 character :: iUplo 398 real(rdp) :: iAlpha, iBeta 399 400 if (present(uplo)) then 401 iUplo = uplo 402 else 403 iUplo = 'L' 404 end if 405 if (present(alpha)) then 406 iAlpha = alpha 407 else 408 iAlpha = 1.0_rdp 409 end if 410 if (present(beta)) then 411 iBeta = beta 412 else 413 iBeta = 0.0_rdp 414 end if 415 416 @:ASSERT(size(y) == size(x)) 417 @:ASSERT(size(a,dim=1) == size(a,dim=2)) 418 @:ASSERT(size(a,dim=1) == size(x)) 419 @:ASSERT(iUplo == 'u' .or. iUplo == 'U' .or. iUplo == 'l' .or. iUplo == 'L') 420 n = size(y) 421 call dsymv( iUplo, n, iAlpha, a, n, x, 1, iBeta, y, 1 ) 422 423 end subroutine symv_dble 424 425 426 !> complex hermitian matrix*vector product 427 subroutine hemv_cmplx(y,a,x,uplo,alpha,beta) 428 429 !> vector 430 complex(rsp), intent(inout) :: y(:) 431 432 !> symmetric matrix 433 complex(rsp), intent(in) :: a(:,:) 434 435 !> vector 436 complex(rsp), intent(in) :: x(:) 437 438 !> optional upper, 'U', or lower 'L' triangle (defaults to lower) 439 character, intent(in), optional :: uplo 440 441 !> optional scaling factor (defaults to 1) 442 complex(rsp), intent(in), optional :: alpha 443 444 !> optional scaling factor (defaults to 0) 445 complex(rsp), intent(in), optional :: beta 446 447 integer :: n 448 character :: iUplo 449 complex(rsp) :: iAlpha, iBeta 450 451 if (present(uplo)) then 452 iUplo = uplo 453 else 454 iUplo = 'L' 455 end if 456 if (present(alpha)) then 457 iAlpha = alpha 458 else 459 iAlpha = cmplx(1.0,0.0,rsp) 460 end if 461 if (present(beta)) then 462 iBeta = beta 463 else 464 iBeta = cmplx(0.0,0.0,rsp) 465 end if 466 467 @:ASSERT(size(y) == size(x)) 468 @:ASSERT(size(a,dim=1) == size(a,dim=2)) 469 @:ASSERT(size(a,dim=1) == size(x)) 470 @:ASSERT(iUplo == 'u' .or. iUplo == 'U' .or. iUplo == 'l' .or. iUplo == 'L') 471 n = size(y) 472 call chemv( iUplo, n, iAlpha, a, n, x, 1, iBeta, y, 1 ) 473 474 end subroutine hemv_cmplx 475 476 477 !> double complex hermitian matrix*vector product 478 subroutine hemv_dblecmplx(y,a,x,uplo,alpha,beta) 479 480 !> vector 481 complex(rdp), intent(inout) :: y(:) 482 483 !> symmetric matrix 484 complex(rdp), intent(in) :: a(:,:) 485 486 !> vector 487 complex(rdp), intent(in) :: x(:) 488 489 !> optional upper, 'U', or lower 'L' triangle (defaults to lower) 490 character, intent(in), optional :: uplo 491 492 !> optional scaling factor (defaults to 1) 493 complex(rdp), intent(in), optional :: alpha 494 495 !> optional scaling factor (defaults to 0) 496 complex(rdp), intent(in), optional :: beta 497 498 integer :: n 499 character :: iUplo 500 complex(rdp) :: iAlpha, iBeta 501 502 if (present(uplo)) then 503 iUplo = uplo 504 else 505 iUplo = 'L' 506 end if 507 if (present(alpha)) then 508 iAlpha = alpha 509 else 510 iAlpha = cmplx(1.0,0.0,rdp) 511 end if 512 if (present(beta)) then 513 iBeta = beta 514 else 515 iBeta = cmplx(0.0,0.0,rdp) 516 end if 517 518 @:ASSERT(size(y) == size(x)) 519 @:ASSERT(size(a,dim=1) == size(a,dim=2)) 520 @:ASSERT(size(a,dim=1) == size(x)) 521 @:ASSERT(iUplo == 'u' .or. iUplo == 'U' .or. iUplo == 'l' .or. iUplo == 'L') 522 n = size(y) 523 call zhemv( iUplo, n, iAlpha, a, n, x, 1, iBeta, y, 1 ) 524 525 end subroutine hemv_dblecmplx 526 527 528 !> real matrix*vector product 529 subroutine gemv_real(y,a,x,alpha,beta,trans) 530 531 !> vector 532 real(rsp), intent(inout) :: y(:) 533 534 !> matrix 535 real(rsp), intent(in) :: a(:,:) 536 537 !> vector 538 real(rsp), intent(in) :: x(:) 539 540 !> optional scaling factor (defaults to 1) 541 real(rsp), intent(in), optional :: alpha 542 543 !> optional scaling factor (defaults to 0) 544 real(rsp), intent(in), optional :: beta 545 546 !> optional transpose (defaults to 'n'), allowed choices are 'n', 'N', 't', 'T', 'c' and 'C' 547 character, intent(in), optional :: trans 548 549 integer :: n, m 550 character :: iTrans 551 real(rsp) :: iAlpha, iBeta 552 553 if (present(trans)) then 554 iTrans = trans 555 else 556 iTrans = 'n' 557 end if 558 if (present(alpha)) then 559 iAlpha = alpha 560 else 561 iAlpha = 1.0_rsp 562 end if 563 if (present(beta)) then 564 iBeta = beta 565 else 566 iBeta = 0.0_rsp 567 end if 568 569 @:ASSERT(iTrans == 'n' .or. iTrans == 'N' .or. iTrans == 't' .or. iTrans == 'T' .or.& 570 & iTrans == 'c' .or. iTrans == 'C') 571 @:ASSERT(((size(a,dim=1) == size(y)) .and. (iTrans == 'n' .or. iTrans == 'N')) .or.& 572 & (size(a,dim=1) == size(x))) 573 @:ASSERT(((size(a,dim=2) == size(x)) .and. (iTrans == 'n' .or. iTrans == 'N')) .or.& 574 & (size(a,dim=2) == size(y))) 575 576 m = size(a,dim=1) 577 n = size(a,dim=2) 578 579 call sgemv( iTrans, m, n, iAlpha, a, m, x, 1, iBeta, y, 1 ) 580 581 end subroutine gemv_real 582 583 584 !> double precision matrix*vector product 585 subroutine gemv_dble(y,a,x,alpha,beta,trans) 586 587 !> vector 588 real(rdp), intent(inout) :: y(:) 589 590 !> matrix 591 real(rdp), intent(in) :: a(:,:) 592 593 !> vector 594 real(rdp), intent(in) :: x(:) 595 596 !> optional scaling factor (defaults to 1) 597 real(rdp), intent(in), optional :: alpha 598 599 !> optional scaling factor (defaults to 0) 600 real(rdp), intent(in), optional :: beta 601 602 !> optional transpose (defaults to 'n'), allowed choices are 'n', 'N', 't', 'T', 'c' and 'C' 603 character, intent(in), optional :: trans 604 605 integer :: n, m 606 character :: iTrans 607 real(rdp) :: iAlpha, iBeta 608 609 if (present(trans)) then 610 iTrans = trans 611 else 612 iTrans = 'n' 613 end if 614 if (present(alpha)) then 615 iAlpha = alpha 616 else 617 iAlpha = 1.0_rdp 618 end if 619 if (present(beta)) then 620 iBeta = beta 621 else 622 iBeta = 0.0_rdp 623 end if 624 625 @:ASSERT(iTrans == 'n' .or. iTrans == 'N' .or. iTrans == 't' .or. iTrans == 'T' .or.& 626 & iTrans == 'c' .or. iTrans == 'C') 627 @:ASSERT(((size(a,dim=1) == size(y)) .and. (iTrans == 'n' .or. iTrans == 'N')) .or.& 628 & (size(a,dim=1) == size(x))) 629 @:ASSERT(((size(a,dim=2) == size(x)) .and. (iTrans == 'n' .or. iTrans == 'N')) .or.& 630 & (size(a,dim=2) == size(y))) 631 632 m = size(a,dim=1) 633 n = size(a,dim=2) 634 635 call dgemv( iTrans, m, n, iAlpha, a, m, x, 1, iBeta, y, 1 ) 636 637 end subroutine gemv_dble 638 639 640 !> real symmetric banded matrix*vector product 641 subroutine sbmv_real(y,ba,x,uplo,alpha,beta) 642 643 !> vector 644 real(rsp), intent(inout) :: y(:) 645 646 !> banded symmetric matrix 647 real(rsp), intent(in) :: ba(:,:) 648 649 !> vector 650 real(rsp), intent(in) :: x(:) 651 652 !> optional upper, 'U', or lower 'L' triangle (defaults to lower) 653 character, intent(in), optional :: uplo 654 655 !> optional scaling factor (defaults to 1) 656 real(rsp), intent(in), optional :: alpha 657 658 !> optional scaling factor (defaults to 0) 659 real(rsp), intent(in), optional :: beta 660 661 integer :: n 662 integer :: k 663 character :: iUplo 664 real(rsp) :: iAlpha, iBeta 665 666 k = size(ba,dim=1) - 1 667 668 if (present(uplo)) then 669 iUplo = uplo 670 else 671 iUplo = 'L' 672 end if 673 if (present(alpha)) then 674 iAlpha = alpha 675 else 676 iAlpha = 1.0_rsp 677 end if 678 if (present(beta)) then 679 iBeta = beta 680 else 681 iBeta = 0.0_rsp 682 end if 683 684 @:ASSERT(size(y) == size(x)) 685 @:ASSERT(size(ba,dim=1) > 0) 686 @:ASSERT(size(ba,dim=2) == size(x)) 687 @:ASSERT(iUplo == 'u' .or. iUplo == 'U' .or. iUplo == 'l' .or. iUplo == 'L') 688 n = size(y) 689 call ssbmv( iUplo, n, k, iAlpha, ba, k+1, x, 1, iBeta, y, 1 ) 690 end subroutine sbmv_real 691 692 693 !> double precision symmetric banded matrix*vector product 694 subroutine sbmv_dble(y,ba,x,uplo,alpha,beta) 695 696 !> vector 697 real(rdp), intent(inout) :: y(:) 698 699 !> banded symmetric matrix 700 real(rdp), intent(in) :: ba(:,:) 701 702 !> vector 703 real(rdp), intent(in) :: x(:) 704 705 !> optional upper, 'U', or lower 'L' triangle (defaults to lower) 706 character, intent(in), optional :: uplo 707 708 !> optional scaling factor (defaults to 1) 709 real(rdp), intent(in), optional :: alpha 710 711 !> optional scaling factor (defaults to 0) 712 real(rdp), intent(in), optional :: beta 713 714 integer :: n 715 integer :: k 716 character :: iUplo 717 real(rdp) :: iAlpha, iBeta 718 719 k = size(ba,dim=1) - 1 720 721 if (present(uplo)) then 722 iUplo = uplo 723 else 724 iUplo = 'L' 725 end if 726 if (present(alpha)) then 727 iAlpha = alpha 728 else 729 iAlpha = 1.0_rdp 730 end if 731 if (present(beta)) then 732 iBeta = beta 733 else 734 iBeta = 0.0_rdp 735 end if 736 737 @:ASSERT(size(y) == size(x)) 738 @:ASSERT(size(ba,dim=1) > 0) 739 @:ASSERT(size(ba,dim=2) == size(x)) 740 @:ASSERT(iUplo == 'u' .or. iUplo == 'U' .or. iUplo == 'l' .or. iUplo == 'L') 741 n = size(y) 742 call dsbmv( iUplo, n, k, iAlpha, ba, k+1, x, 1, iBeta, y, 1 ) 743 end subroutine sbmv_dble 744 745 746 !> real precision symmetric matrix * general matrix multiply 747 subroutine symm_real(C,side,A,B,uplo,alpha,beta,m,n) 748 749 !> general matrix output 750 real(rsp), intent(inout) :: C(:,:) 751 752 !> symmetric matrix on 'l'eft or 'r'ight , where A is symmetric and B is general SIDE = 'L' or 753 !> 'l' C := alpha*A*B + beta*C, SIDE = 'R' or 'r' C := alpha*B*A + beta*C 754 character, intent(in) :: side 755 756 !> symmetric matrix, size 757 real(rsp), intent(in) :: A(:,:) 758 759 !> general matrix 760 real(rsp), intent(in) :: B(:,:) 761 762 !> is an 'U'pper or 'L'ower triangle matrix, defaults to lower 763 character, intent(in), optional :: uplo 764 765 !> defaults to 1 if not set 766 real(rsp), intent(in), optional :: alpha 767 768 !> defaults to 0 if not set 769 real(rsp), intent(in), optional :: beta 770 771 !> specifies the number of rows of the matrix C 772 integer, intent(in), optional :: m 773 774 !> specifies the number of columns of the matrix C 775 integer, intent(in), optional :: n 776 777 integer :: lda, ldb, ldc, ka, im, in 778 character :: iUplo 779 real(rsp) :: iAlpha, iBeta 780 781 if (present(uplo)) then 782 iUplo = uplo 783 else 784 iUplo = 'L' 785 end if 786 if (present(alpha)) then 787 iAlpha = alpha 788 else 789 iAlpha = 1.0_rsp 790 end if 791 if (present(beta)) then 792 iBeta = beta 793 else 794 iBeta = 0.0_rsp 795 end if 796 797 lda = size(a,dim=1) 798 ldb = size(b,dim=1) 799 ldc = size(c,dim=1) 800 801 @:ASSERT(iUplo == 'u' .or. iUplo == 'U' .or. iUplo == 'l' .or. iUplo == 'L') 802 @:ASSERT(side == 'r' .or. side == 'R' .or. side == 'l' .or. side == 'L') 803 804 if (present(n)) then 805 in = n 806 else 807 in = size(C,dim=2) 808 end if 809 if (present(m)) then 810 im = m 811 else 812 im = size(C,dim=1) 813 end if 814 if (iUplo=='l'.or.iUplo=='L') then 815 ka = im 816 else 817 ka = in 818 end if 819 820 @:ASSERT(im>0) 821 @:ASSERT(in>0) 822 @:ASSERT((lda>=im).and.(iUplo=='l'.or.iUplo=='L').or.(lda>=in)) 823 @:ASSERT(size(a,dim=2)>=ka) 824 @:ASSERT(ldb>=im) 825 @:ASSERT(ldc>=im) 826 @:ASSERT(size(B,dim=2)>=in) 827 @:ASSERT(size(C,dim=2)>=in) 828 829 call ssymm ( side, iUplo, im, in, iAlpha, A, lda, B, ldb, iBeta, C, ldc ) 830 831 end subroutine symm_real 832 833 834 !> double precision symmetric matrix * general matrix multiply 835 subroutine symm_dble(C,side,A,B,uplo,alpha,beta,m,n) 836 837 !> general matrix output 838 real(rdp), intent(inout) :: C(:,:) 839 840 !> symmetric matrix on 'l'eft or 'r'ight , where A is symmetric and B is general SIDE = 'L' or 841 !> 'l' C := alpha*A*B + beta*C, SIDE = 'R' or 'r' C := alpha*B*A + beta*C 842 character, intent(in) :: side 843 844 !> symmetric matrix, size 845 real(rdp), intent(in) :: A(:,:) 846 847 !> general matrix 848 real(rdp), intent(in) :: B(:,:) 849 850 !> is an 'U'pper or 'L'ower triangle matrix, defaults to lower 851 character, intent(in), optional :: uplo 852 853 !> defaults to 1 if not set 854 real(rdp), intent(in), optional :: alpha 855 856 !> defaults to 0 if not set 857 real(rdp), intent(in), optional :: beta 858 859 !> specifies the number of rows of the matrix C 860 integer, intent(in), optional :: m 861 862 !> specifies the number of columns of the matrix C 863 integer, intent(in), optional :: n 864 865 integer :: lda, ldb, ldc, ka, im, in 866 character :: iUplo 867 real(rdp) :: iAlpha, iBeta 868 869 if (present(uplo)) then 870 iUplo = uplo 871 else 872 iUplo = 'L' 873 end if 874 if (present(alpha)) then 875 iAlpha = alpha 876 else 877 iAlpha = 1.0_rdp 878 end if 879 if (present(beta)) then 880 iBeta = beta 881 else 882 iBeta = 0.0_rdp 883 end if 884 885 lda = size(a,dim=1) 886 ldb = size(b,dim=1) 887 ldc = size(c,dim=1) 888 889 @:ASSERT(iUplo == 'u' .or. iUplo == 'U' .or. iUplo == 'l' .or. iUplo == 'L') 890 @:ASSERT(side == 'r' .or. side == 'R' .or. side == 'l' .or. side == 'L') 891 892 if (present(n)) then 893 in = n 894 else 895 in = size(C,dim=2) 896 end if 897 if (present(m)) then 898 im = m 899 else 900 im = size(C,dim=1) 901 end if 902 if (iUplo=='l'.or.iUplo=='L') then 903 ka = im 904 else 905 ka = in 906 end if 907 908 @:ASSERT(im>0) 909 @:ASSERT(in>0) 910 @:ASSERT((lda>=im).and.(iUplo=='l'.or.iUplo=='L').or.(lda>=in)) 911 @:ASSERT(size(a,dim=2)>=ka) 912 @:ASSERT(ldb>=im) 913 @:ASSERT(ldc>=im) 914 @:ASSERT(size(B,dim=2)>=in) 915 @:ASSERT(size(C,dim=2)>=in) 916 917 call dsymm ( side, iUplo, im, in, iAlpha, A, lda, B, ldb, iBeta, C, ldc ) 918 919 end subroutine symm_dble 920 921 922 !> real matrix*matrix product 923 subroutine gemm_real(C,A,B,alpha,beta,transA,transB,n,m,k) 924 925 !> general matrix output 926 real(rsp), intent(inout) :: C(:,:) 927 928 !> symmetric matrix 929 real(rsp), intent(in) :: A(:,:) 930 931 !> general matrix 932 real(rsp), intent(in) :: B(:,:) 933 934 !> defaults to 1 if not set 935 real(rsp), intent(in), optional :: alpha 936 937 !> defaults to 0 if not set 938 real(rsp), intent(in), optional :: beta 939 940 !> optional transpose of A matrix (defaults to 'n'), allowed choices are 'n', 'N', 't', 'T', 'c' 941 !> and 'C' 942 character, intent(in), optional :: transA 943 944 !> optional transpose of B matrix (defaults to 'n'), allowed choices are 'n', 'N', 't', 'T', 'c' 945 !> and 'C' 946 character, intent(in), optional :: transB 947 948 !> specifies the number of columns of the matrix C 949 integer, intent(in), optional :: n 950 951 !> specifies the number of rows of the matrix C 952 integer, intent(in), optional :: m 953 954 !> specifies the internal number of elements in Op(A)_ik Op(B)_kj 955 integer, intent(in), optional :: k 956 957 integer :: lda, ldb, ldc 958 integer :: in, im, ik 959 character :: iTransA, iTransB 960 real(rsp) :: iAlpha, iBeta 961 962 if (present(transA)) then 963 iTransA = transA 964 else 965 iTransA = 'n' 966 end if 967 if (present(transB)) then 968 iTransB = transB 969 else 970 iTransB = 'n' 971 end if 972 973 @:ASSERT(iTransA == 'n' .or. iTransA == 'N' .or. iTransA == 't'& 974 & .or. iTransA == 'T' .or. iTransA == 'c' .or. iTransA == 'C') 975 @:ASSERT(iTransB == 'n' .or. iTransB == 'N' .or. iTransB == 't'& 976 & .or. iTransB == 'T' .or. iTransB == 'c' .or. iTransB == 'C') 977 978 if (present(alpha)) then 979 iAlpha = alpha 980 else 981 iAlpha = 1.0_rsp 982 end if 983 if (present(beta)) then 984 iBeta = beta 985 else 986 iBeta = 0.0_rsp 987 end if 988 989 lda = size(a,dim=1) 990 ldb = size(b,dim=1) 991 ldc = size(c,dim=1) 992 993 if (present(m)) then 994 im = m 995 else 996 if (iTransA == 'n' .or. iTransA == 'N') then 997 im = size(A,dim=1) 998 else 999 im = size(A,dim=2) 1000 end if 1001 end if 1002 if (present(n)) then 1003 in = n 1004 else 1005 in = size(c,dim=2) 1006 end if 1007 if (present(k)) then 1008 ik = k 1009 else 1010 if (iTransA == 'n' .or. iTransA == 'N') then 1011 ik = size(A,dim=2) 1012 else 1013 ik = size(A,dim=1) 1014 end if 1015 end if 1016 1017 @:ASSERT(im>0) 1018 @:ASSERT(in>0) 1019 @:ASSERT(ik>0) 1020 @:ASSERT(((lda>=im).and.(iTransA == 'n' .or. iTransA == 'N'))& 1021 & .or. (size(a,dim=2)>=im)) 1022 @:ASSERT(ldc>=im) 1023 @:ASSERT(((size(b,dim=2)>=in).and.(iTransB == 'n' .or. iTransB == 'N'))& 1024 & .or. (ldb>=in)) 1025 @:ASSERT(size(c,dim=2)>=in) 1026 @:ASSERT(((size(a,dim=2)>=ik).and.(iTransA == 'n' .or. iTransA == 'N'))& 1027 & .or. (lda>=ik)) 1028 @:ASSERT(((ldb>=ik).and.(iTransB == 'n' .or. iTransB == 'N'))& 1029 & .or. (size(b,dim=2)>=ik)) 1030 1031 call sgemm(iTransA,iTransB,im,in,ik,iAlpha,A,lda,B,ldb,iBeta,C,ldc) 1032 1033 end subroutine gemm_real 1034 1035 1036 !> Double precision matrix*matrix product 1037 subroutine gemm_dble(C,A,B,alpha,beta,transA,transB,n,m,k) 1038 1039 !> general matrix output 1040 real(rdp), intent(inout) :: C(:,:) 1041 1042 !> symmetric matrix 1043 real(rdp), intent(in) :: A(:,:) 1044 1045 !> general matrix 1046 real(rdp), intent(in) :: B(:,:) 1047 1048 !> defaults to 1 if not set 1049 real(rdp), intent(in), optional :: alpha 1050 1051 !> defaults to 0 if not set 1052 real(rdp), intent(in), optional :: beta 1053 1054 !> optional transpose of A matrix (defaults to 'n'), allowed choices are 'n', 'N', 't', 'T', 'c' 1055 !> and 'C' 1056 character, intent(in), optional :: transA 1057 1058 !> optional transpose of B matrix (defaults to 'n'), allowed choices are 'n', 'N', 't', 'T', 'c' 1059 !> and 'C' 1060 character, intent(in), optional :: transB 1061 1062 !> specifies the number of columns of the matrix C 1063 integer, intent(in), optional :: n 1064 1065 !> specifies the number of rows of the matrix C 1066 integer, intent(in), optional :: m 1067 1068 !> specifies the internal number of elements in Op(A)_ik Op(B)_kj 1069 integer, intent(in), optional :: k 1070 1071 integer :: lda, ldb, ldc 1072 integer :: in, im, ik 1073 character :: iTransA, iTransB 1074 real(rdp) :: iAlpha, iBeta 1075 1076 if (present(transA)) then 1077 iTransA = transA 1078 else 1079 iTransA = 'n' 1080 end if 1081 if (present(transB)) then 1082 iTransB = transB 1083 else 1084 iTransB = 'n' 1085 end if 1086 1087 @:ASSERT(iTransA == 'n' .or. iTransA == 'N' .or. iTransA == 't'& 1088 & .or. iTransA == 'T' .or. iTransA == 'c' .or. iTransA == 'C') 1089 @:ASSERT(iTransB == 'n' .or. iTransB == 'N' .or. iTransB == 't'& 1090 & .or. iTransB == 'T' .or. iTransB == 'c' .or. iTransB == 'C') 1091 1092 if (present(alpha)) then 1093 iAlpha = alpha 1094 else 1095 iAlpha = 1.0_rdp 1096 end if 1097 if (present(beta)) then 1098 iBeta = beta 1099 else 1100 iBeta = 0.0_rdp 1101 end if 1102 1103 lda = size(a,dim=1) 1104 ldb = size(b,dim=1) 1105 ldc = size(c,dim=1) 1106 1107 if (present(m)) then 1108 im = m 1109 else 1110 if (iTransA == 'n' .or. iTransA == 'N') then 1111 im = size(A,dim=1) 1112 else 1113 im = size(A,dim=2) 1114 end if 1115 end if 1116 if (present(n)) then 1117 in = n 1118 else 1119 in = size(c,dim=2) 1120 end if 1121 if (present(k)) then 1122 ik = k 1123 else 1124 if (iTransA == 'n' .or. iTransA == 'N') then 1125 ik = size(A,dim=2) 1126 else 1127 ik = size(A,dim=1) 1128 end if 1129 end if 1130 1131 @:ASSERT(im>0) 1132 @:ASSERT(in>0) 1133 @:ASSERT(ik>0) 1134 @:ASSERT(((lda>=im).and.(iTransA == 'n' .or. iTransA == 'N'))& 1135 & .or. (size(a,dim=2)>=im)) 1136 @:ASSERT(ldc>=im) 1137 @:ASSERT(((size(b,dim=2)>=in).and.(iTransB == 'n' .or. iTransB == 'N'))& 1138 & .or. (ldb>=in)) 1139 @:ASSERT(size(c,dim=2)>=in) 1140 @:ASSERT(((size(a,dim=2)>=ik).and.(iTransA == 'n' .or. iTransA == 'N'))& 1141 & .or. (lda>=ik)) 1142 @:ASSERT(((ldb>=ik).and.(iTransB == 'n' .or. iTransB == 'N'))& 1143 & .or. (size(b,dim=2)>=ik)) 1144 1145 call dgemm(iTransA,iTransB,im,in,ik,iAlpha,A,lda,B,ldb,iBeta,C,ldc) 1146 1147 end subroutine gemm_dble 1148 1149 1150 !> complex matrix*matrix product 1151 subroutine gemm_cmplx(C,A,B,alpha,beta,transA,transB,n,m,k) 1152 1153 !> general matrix output 1154 complex(rsp), intent(inout) :: C(:,:) 1155 1156 !> symmetric matrix 1157 complex(rsp), intent(in) :: A(:,:) 1158 1159 !> general matrix 1160 complex(rsp), intent(in) :: B(:,:) 1161 1162 !> defaults to 1 if not set 1163 complex(rsp), intent(in), optional :: alpha 1164 1165 !> defaults to 0 if not set 1166 complex(rsp), intent(in), optional :: beta 1167 1168 !> optional transpose of A matrix (defaults to 'n'), allowed choices are 'n', 'N', 't', 'T', 'c' 1169 !> and 'C' 1170 character, intent(in), optional :: transA 1171 1172 !> optional transpose of B matrix (defaults to 'n'), allowed choices are 'n', 'N', 't', 'T', 'c' 1173 !> and 'C' 1174 character, intent(in), optional :: transB 1175 1176 !> specifies the number of columns of the matrix C 1177 integer, intent(in), optional :: n 1178 1179 !> specifies the number of rows of the matrix C 1180 integer, intent(in), optional :: m 1181 1182 !> specifies the internal number of elements in Op(A)_ik Op(B)_kj 1183 integer, intent(in), optional :: k 1184 1185 integer :: lda, ldb, ldc 1186 integer :: in, im, ik 1187 character :: iTransA, iTransB 1188 complex(rsp) :: iAlpha, iBeta 1189 1190 if (present(transA)) then 1191 iTransA = transA 1192 else 1193 iTransA = 'n' 1194 end if 1195 if (present(transB)) then 1196 iTransB = transB 1197 else 1198 iTransB = 'n' 1199 end if 1200 1201 @:ASSERT(iTransA == 'n' .or. iTransA == 'N' .or. iTransA == 't'& 1202 & .or. iTransA == 'T' .or. iTransA == 'c' .or. iTransA == 'C') 1203 @:ASSERT(iTransB == 'n' .or. iTransB == 'N' .or. iTransB == 't'& 1204 & .or. iTransB == 'T' .or. iTransB == 'c' .or. iTransB == 'C') 1205 1206 if (present(alpha)) then 1207 iAlpha = alpha 1208 else 1209 iAlpha = cmplx(1.0,0.0,rsp) 1210 end if 1211 if (present(beta)) then 1212 iBeta = beta 1213 else 1214 iBeta = cmplx(0.0,0.0,rsp) 1215 end if 1216 1217 lda = size(a,dim=1) 1218 ldb = size(b,dim=1) 1219 ldc = size(c,dim=1) 1220 1221 if (present(m)) then 1222 im = m 1223 else 1224 if (iTransA == 'n' .or. iTransA == 'N') then 1225 im = size(A,dim=1) 1226 else 1227 im = size(A,dim=2) 1228 end if 1229 end if 1230 if (present(n)) then 1231 in = n 1232 else 1233 in = size(c,dim=2) 1234 end if 1235 if (present(k)) then 1236 ik = k 1237 else 1238 if (iTransA == 'n' .or. iTransA == 'N') then 1239 ik = size(A,dim=2) 1240 else 1241 ik = size(A,dim=1) 1242 end if 1243 end if 1244 1245 @:ASSERT(im>0) 1246 @:ASSERT(in>0) 1247 @:ASSERT(ik>0) 1248 @:ASSERT(((lda>=im).and.(iTransA == 'n' .or. iTransA == 'N'))& 1249 & .or. (size(a,dim=2)>=im)) 1250 @:ASSERT(ldc>=im) 1251 @:ASSERT(((size(b,dim=2)>=in).and.(iTransB == 'n' .or. iTransB == 'N'))& 1252 & .or. (ldb>=in)) 1253 @:ASSERT(size(c,dim=2)>=in) 1254 @:ASSERT(((size(a,dim=2)>=ik).and.(iTransA == 'n' .or. iTransA == 'N'))& 1255 & .or. (lda>=ik)) 1256 @:ASSERT(((ldb>=ik).and.(iTransB == 'n' .or. iTransB == 'N'))& 1257 & .or. (size(b,dim=2)>=ik)) 1258 1259 call cgemm(iTransA,iTransB,im,in,ik,iAlpha,A,lda,B,ldb,iBeta,C,ldc) 1260 1261 end subroutine gemm_cmplx 1262 1263 1264 !> Double precision matrix*matrix product 1265 subroutine gemm_dblecmplx(C,A,B,alpha,beta,transA,transB,n,m,k) 1266 1267 !> general matrix output 1268 complex(rdp), intent(inout) :: C(:,:) 1269 1270 !> symmetric matrix 1271 complex(rdp), intent(in) :: A(:,:) 1272 1273 !> general matrix 1274 complex(rdp), intent(in) :: B(:,:) 1275 1276 !> defaults to 1 if not set 1277 complex(rdp), intent(in), optional :: alpha 1278 1279 !> defaults to 0 if not set 1280 complex(rdp), intent(in), optional :: beta 1281 1282 !> optional transpose of A matrix (defaults to 'n'), allowed choices are 'n', 'N', 't', 'T', 'c' 1283 !> and 'C' 1284 character, intent(in), optional :: transA 1285 1286 !> optional transpose of B matrix (defaults to 'n'), allowed choices are 'n', 'N', 't', 'T', 'c' 1287 !> and 'C' 1288 character, intent(in), optional :: transB 1289 1290 !> specifies the number of columns of the matrix C 1291 integer, intent(in), optional :: n 1292 1293 !> specifies the number of rows of the matrix C 1294 integer, intent(in), optional :: m 1295 1296 !> specifies the internal number of elements in Op(A)_ik Op(B)_kj 1297 integer, intent(in), optional :: k 1298 1299 integer :: lda, ldb, ldc 1300 integer :: in, im, ik 1301 character :: iTransA, iTransB 1302 complex(rdp) :: iAlpha, iBeta 1303 1304 if (present(transA)) then 1305 iTransA = transA 1306 else 1307 iTransA = 'n' 1308 end if 1309 if (present(transB)) then 1310 iTransB = transB 1311 else 1312 iTransB = 'n' 1313 end if 1314 1315 @:ASSERT(iTransA == 'n' .or. iTransA == 'N' .or. iTransA == 't'& 1316 & .or. iTransA == 'T' .or. iTransA == 'c' .or. iTransA == 'C') 1317 @:ASSERT(iTransB == 'n' .or. iTransB == 'N' .or. iTransB == 't'& 1318 & .or. iTransB == 'T' .or. iTransB == 'c' .or. iTransB == 'C') 1319 1320 if (present(alpha)) then 1321 iAlpha = alpha 1322 else 1323 iAlpha = cmplx(1.0,0.0,rdp) 1324 end if 1325 if (present(beta)) then 1326 iBeta = beta 1327 else 1328 iBeta = cmplx(0.0,0.0,rdp) 1329 end if 1330 1331 lda = size(a,dim=1) 1332 ldb = size(b,dim=1) 1333 ldc = size(c,dim=1) 1334 1335 if (present(m)) then 1336 im = m 1337 else 1338 if (iTransA == 'n' .or. iTransA == 'N') then 1339 im = size(A,dim=1) 1340 else 1341 im = size(A,dim=2) 1342 end if 1343 end if 1344 if (present(n)) then 1345 in = n 1346 else 1347 in = size(c,dim=2) 1348 end if 1349 if (present(k)) then 1350 ik = k 1351 else 1352 if (iTransA == 'n' .or. iTransA == 'N') then 1353 ik = size(A,dim=2) 1354 else 1355 ik = size(A,dim=1) 1356 end if 1357 end if 1358 1359 @:ASSERT(im>0) 1360 @:ASSERT(in>0) 1361 @:ASSERT(ik>0) 1362 @:ASSERT(((lda>=im).and.(iTransA == 'n' .or. iTransA == 'N'))& 1363 & .or. (size(a,dim=2)>=im)) 1364 @:ASSERT(ldc>=im) 1365 @:ASSERT(((size(b,dim=2)>=in).and.(iTransB == 'n' .or. iTransB == 'N'))& 1366 & .or. (ldb>=in)) 1367 @:ASSERT(size(c,dim=2)>=in) 1368 @:ASSERT(((size(a,dim=2)>=ik).and.(iTransA == 'n' .or. iTransA == 'N'))& 1369 & .or. (lda>=ik)) 1370 @:ASSERT(((ldb>=ik).and.(iTransB == 'n' .or. iTransB == 'N'))& 1371 & .or. (size(b,dim=2)>=ik)) 1372 1373 call zgemm(iTransA,iTransB,im,in,ik,iAlpha,A,lda,B,ldb,iBeta,C,ldc) 1374 1375 end subroutine gemm_dblecmplx 1376 1377 1378 !> real rank-k update 1379 subroutine herk_real(C,A,alpha,beta,uplo,trans,n,k) 1380 1381 !> contains the matrix to be updated 1382 real(rsp), intent(inout) :: C(:,:) 1383 1384 !> contains the matrix to update 1385 real(rsp), intent(in) :: A(:,:) 1386 1387 !> scaling value for the update contribution, defaults to 1 1388 real(rsp), intent(in), optional :: alpha 1389 1390 !> scaling value for the original C, defaults to 0 1391 real(rsp), intent(in), optional :: beta 1392 1393 !> optional upper, 'U', or lower 'L' triangle, defaults to lower 1394 character, intent(in), optional :: uplo 1395 1396 !> optional transpose (defaults to 'n'), allowed choices are 'n', 'N', 't', 'T' (and 'C' or 'c' 1397 !> for the real cases) 1398 character, intent(in), optional :: trans 1399 1400 !> order of the matrix C 1401 integer, intent(in), optional :: n 1402 1403 !> internal order of A summation 1404 integer, intent(in), optional :: k 1405 1406 integer :: lda, ldc 1407 integer :: in, ik 1408 character :: iTrans, iUplo 1409 real(rsp) :: iAlpha, iBeta 1410 1411 if (present(uplo)) then 1412 iUplo = uplo 1413 else 1414 iUplo = 'L' 1415 end if 1416 @:ASSERT(iUplo == 'u' .or. iUplo == 'U' .or. iUplo == 'l' .or. iUplo == 'L') 1417 1418 if (present(trans)) then 1419 iTrans = trans 1420 else 1421 iTrans = 'n' 1422 end if 1423 1424 @:ASSERT(iTrans == 'n' .or. iTrans == 'N' .or. iTrans == 't'& 1425 & .or. iTrans == 'T' .or. iTrans == 'c' .or. iTrans == 'C') 1426 1427 if (present(alpha)) then 1428 iAlpha = alpha 1429 else 1430 iAlpha = 1.0_rsp 1431 end if 1432 if (present(beta)) then 1433 iBeta = beta 1434 else 1435 iBeta = 0.0_rsp 1436 end if 1437 1438 lda = size(a,dim=1) 1439 ldc = size(c,dim=1) 1440 1441 if (present(n)) then 1442 in = n 1443 else 1444 in = size(c,dim=2) 1445 end if 1446 if (present(k)) then 1447 ik = k 1448 else 1449 if (iTrans == 'n' .or. iTrans == 'N') then 1450 ik = size(A,dim=2) 1451 else 1452 ik = size(A,dim=1) 1453 end if 1454 end if 1455 1456 @:ASSERT(in>0) 1457 @:ASSERT(ik>0) 1458 @:ASSERT(((size(a,dim=2)>=in).and.(iTrans == 'n' .or. iTrans == 'N'))& 1459 & .or. (lda>=in)) 1460 @:ASSERT(size(c,dim=2)>=in) 1461 @:ASSERT(((size(a,dim=2)>=ik).and.(iTrans == 'n' .or. iTrans == 'N'))& 1462 & .or. (lda>=ik)) 1463 1464 call ssyrk(iUplo, iTrans, in, ik, iAlpha, A, lda, iBeta, C, ldc ) 1465 1466 end subroutine herk_real 1467 1468 1469 !> Double precision rank-k update 1470 subroutine herk_dble(C,A,alpha,beta,uplo,trans,n,k) 1471 1472 !> contains the matrix to be updated 1473 real(rdp), intent(inout) :: C(:,:) 1474 1475 !> contains the matrix to update 1476 real(rdp), intent(in) :: A(:,:) 1477 1478 !> scaling value for the update contribution, defaults to 1 1479 real(rdp), intent(in), optional :: alpha 1480 1481 !> scaling value for the original C, defaults to 0 1482 real(rdp), intent(in), optional :: beta 1483 1484 !> optional upper, 'U', or lower 'L' triangle, defaults to lower 1485 character, intent(in), optional :: uplo 1486 1487 !> optional transpose (defaults to 'n'), allowed choices are 'n', 'N', 't', 'T' (and 'C' or 'c' 1488 !> for the real cases) 1489 character, intent(in), optional :: trans 1490 1491 !> order of the matrix C 1492 integer, intent(in), optional :: n 1493 1494 !> internal order of A summation 1495 integer, intent(in), optional :: k 1496 1497 integer :: lda, ldc 1498 integer :: in, ik 1499 character :: iTrans, iUplo 1500 real(rdp) :: iAlpha, iBeta 1501 1502 if (present(uplo)) then 1503 iUplo = uplo 1504 else 1505 iUplo = 'L' 1506 end if 1507 @:ASSERT(iUplo == 'u' .or. iUplo == 'U' .or. iUplo == 'l' .or. iUplo == 'L') 1508 1509 if (present(trans)) then 1510 iTrans = trans 1511 else 1512 iTrans = 'n' 1513 end if 1514 1515 @:ASSERT(iTrans == 'n' .or. iTrans == 'N' .or. iTrans == 't'& 1516 & .or. iTrans == 'T' .or. iTrans == 'c' .or. iTrans == 'C') 1517 1518 if (present(alpha)) then 1519 iAlpha = alpha 1520 else 1521 iAlpha = 1.0_rdp 1522 end if 1523 if (present(beta)) then 1524 iBeta = beta 1525 else 1526 iBeta = 0.0_rdp 1527 end if 1528 1529 lda = size(a,dim=1) 1530 ldc = size(c,dim=1) 1531 1532 if (present(n)) then 1533 in = n 1534 else 1535 in = size(c,dim=2) 1536 end if 1537 if (present(k)) then 1538 ik = k 1539 else 1540 if (iTrans == 'n' .or. iTrans == 'N') then 1541 ik = size(A,dim=2) 1542 else 1543 ik = size(A,dim=1) 1544 end if 1545 end if 1546 1547 @:ASSERT(in>0) 1548 @:ASSERT(ik>0) 1549 @:ASSERT(((size(a,dim=2)>=in).and.(iTrans == 'n' .or. iTrans == 'N'))& 1550 & .or. (lda>=in)) 1551 @:ASSERT(size(c,dim=2)>=in) 1552 @:ASSERT(((size(a,dim=2)>=ik).and.(iTrans == 'n' .or. iTrans == 'N'))& 1553 & .or. (lda>=ik)) 1554 1555 call dsyrk(iUplo, iTrans, in, ik, iAlpha, A, lda, iBeta, C, ldc ) 1556 1557 end subroutine herk_dble 1558 1559 1560 !> complex rank-k update 1561 subroutine herk_cmplx(C,A,alpha,beta,uplo,trans,n,k) 1562 1563 !> contains the matrix to be updated 1564 complex(rsp), intent(inout) :: C(:,:) 1565 1566 !> contains the matrix to update 1567 complex(rsp), intent(in) :: A(:,:) 1568 1569 !> scaling value for the update contribution, defaults to 1 1570 real(rsp), intent(in), optional :: alpha 1571 1572 !> scaling value for the original C, defaults to 0 1573 real(rsp), intent(in), optional :: beta 1574 1575 !> optional upper, 'U', or lower 'L' triangle, defaults to lower 1576 character, intent(in), optional :: uplo 1577 1578 !> optional transpose (defaults to 'n'), allowed choices are 'n', 'N', 't', 'T' (and 'C' or 'c' 1579 !> for the real cases) 1580 character, intent(in), optional :: trans 1581 1582 !> order of the matrix C 1583 integer, intent(in), optional :: n 1584 1585 !> internal order of A summation 1586 integer, intent(in), optional :: k 1587 1588 integer :: lda, ldc 1589 integer :: in, ik 1590 character :: iTrans, iUplo 1591 real(rsp) :: iAlpha, iBeta 1592 1593 if (present(uplo)) then 1594 iUplo = uplo 1595 else 1596 iUplo = 'L' 1597 end if 1598 @:ASSERT(iUplo == 'u' .or. iUplo == 'U' .or. iUplo == 'l' .or. iUplo == 'L') 1599 1600 if (present(trans)) then 1601 iTrans = trans 1602 else 1603 iTrans = 'n' 1604 end if 1605 1606 @:ASSERT(iTrans == 'n' .or. iTrans == 'N' .or. iTrans == 't' .or. iTrans == 'T') 1607 1608 if (present(alpha)) then 1609 iAlpha = alpha 1610 else 1611 iAlpha = 1.0_rsp 1612 end if 1613 if (present(beta)) then 1614 iBeta = beta 1615 else 1616 iBeta = 0.0_rsp 1617 end if 1618 1619 lda = size(a,dim=1) 1620 ldc = size(c,dim=1) 1621 1622 if (present(n)) then 1623 in = n 1624 else 1625 in = size(c,dim=2) 1626 end if 1627 if (present(k)) then 1628 ik = k 1629 else 1630 if (iTrans == 'n' .or. iTrans == 'N') then 1631 ik = size(A,dim=2) 1632 else 1633 ik = size(A,dim=1) 1634 end if 1635 end if 1636 1637 @:ASSERT(in>0) 1638 @:ASSERT(ik>0) 1639 @:ASSERT(((size(a,dim=2)>=in).and.(iTrans == 'n' .or. iTrans == 'N')) .or. (lda>=in)) 1640 @:ASSERT(size(c,dim=2)>=in) 1641 @:ASSERT(((size(a,dim=2)>=ik).and.(iTrans == 'n' .or. iTrans == 'N')) .or. (lda>=ik)) 1642 1643 call cherk(iUplo, iTrans, in, ik, iAlpha, A, lda, iBeta, C, ldc ) 1644 1645 end subroutine herk_cmplx 1646 1647 1648 !> Double complex rank-k update 1649 subroutine herk_dblecmplx(C,A,alpha,beta,uplo,trans,n,k) 1650 1651 !> contains the matrix to be updated 1652 complex(rdp), intent(inout) :: C(:,:) 1653 1654 !> contains the matrix to update 1655 complex(rdp), intent(in) :: A(:,:) 1656 1657 !> scaling value for the update contribution, defaults to 1 1658 real(rdp), intent(in), optional :: alpha 1659 1660 !> scaling value for the original C, defaults to 0 1661 real(rdp), intent(in), optional :: beta 1662 1663 !> optional upper, 'U', or lower 'L' triangle, defaults to lower 1664 character, intent(in), optional :: uplo 1665 1666 !> optional transpose (defaults to 'n'), allowed choices are 'n', 'N', 't', 'T' (and 'C' or 'c' 1667 !> for the real cases) 1668 character, intent(in), optional :: trans 1669 1670 !> order of the matrix C 1671 integer, intent(in), optional :: n 1672 1673 !> internal order of A summation 1674 integer, intent(in), optional :: k 1675 1676 integer :: lda, ldc 1677 integer :: in, ik 1678 character :: iTrans, iUplo 1679 real(rdp) :: iAlpha, iBeta 1680 1681 if (present(uplo)) then 1682 iUplo = uplo 1683 else 1684 iUplo = 'L' 1685 end if 1686 @:ASSERT(iUplo == 'u' .or. iUplo == 'U' .or. iUplo == 'l' .or. iUplo == 'L') 1687 1688 if (present(trans)) then 1689 iTrans = trans 1690 else 1691 iTrans = 'n' 1692 end if 1693 1694 @:ASSERT(iTrans == 'n' .or. iTrans == 'N' .or. iTrans == 't' .or. iTrans == 'T') 1695 1696 if (present(alpha)) then 1697 iAlpha = alpha 1698 else 1699 iAlpha = 1.0_rdp 1700 end if 1701 if (present(beta)) then 1702 iBeta = beta 1703 else 1704 iBeta = 0.0_rdp 1705 end if 1706 1707 lda = size(a,dim=1) 1708 ldc = size(c,dim=1) 1709 1710 if (present(n)) then 1711 in = n 1712 else 1713 in = size(c,dim=2) 1714 end if 1715 if (present(k)) then 1716 ik = k 1717 else 1718 if (iTrans == 'n' .or. iTrans == 'N') then 1719 ik = size(A,dim=2) 1720 else 1721 ik = size(A,dim=1) 1722 end if 1723 end if 1724 1725 @:ASSERT(in>0) 1726 @:ASSERT(ik>0) 1727 @:ASSERT(((size(a,dim=2)>=in).and.(iTrans == 'n' .or. iTrans == 'N')) .or. (lda>=in)) 1728 @:ASSERT(size(c,dim=2)>=in) 1729 @:ASSERT(((size(a,dim=2)>=ik).and.(iTrans == 'n' .or. iTrans == 'N')) .or. (lda>=ik)) 1730 1731 call zherk(iUplo, iTrans, in, ik, iAlpha, A, lda, iBeta, C, ldc ) 1732 1733 end subroutine herk_dblecmplx 1734 1735 1736 !> single precision hermitian matrix * general matrix multiply 1737 subroutine hemm_cmplx(C, side, A, B, uplo, alpha, beta, m, n) 1738 1739 !> general matrix output 1740 complex(rsp), intent(inout) :: C(:,:) 1741 1742 !> symmetric matrix on 'l'eft or 'r'ight , where A is symmetric and B is general SIDE = 'L' or 1743 !> 'l' C := alpha*A*B + beta*C, SIDE = 'R' or 'r' C := alpha*B*A + beta*C 1744 character, intent(in) :: side 1745 1746 !> hermitian matrix 1747 complex(rsp), intent(in) :: A(:,:) 1748 1749 !> general matrix 1750 complex(rsp), intent(in) :: B(:,:) 1751 1752 !> A is an 'U'pper or 'L'ower triangle matrix, defaults to lower 1753 character, intent(in), optional :: uplo 1754 1755 !> defaults to 1 if not set 1756 complex(rsp), intent(in), optional :: alpha 1757 1758 !> defaults to 0 if not set 1759 complex(rsp), intent(in), optional :: beta 1760 1761 !> specifies the number of rows of the matrix C 1762 integer, intent(in), optional :: m 1763 1764 !> specifies the number of columns of the matrix C 1765 integer, intent(in), optional :: n 1766 1767 integer :: lda, ldb, ldc, ka, im, in 1768 character :: iUplo 1769 complex(rsp) :: iAlpha, iBeta 1770 1771 if (present(uplo)) then 1772 iUplo = uplo 1773 else 1774 iUplo = 'L' 1775 end if 1776 if (present(alpha)) then 1777 iAlpha = alpha 1778 else 1779 iAlpha = (1.0_rsp, 0.0_rsp) 1780 end if 1781 if (present(beta)) then 1782 iBeta = beta 1783 else 1784 iBeta = (0.0_rsp, 0.0_rsp) 1785 end if 1786 1787 lda = size(a,dim=1) 1788 ldb = size(b,dim=1) 1789 ldc = size(c,dim=1) 1790 1791 @:ASSERT(iUplo == 'u' .or. iUplo == 'U' .or. iUplo == 'l' .or. iUplo == 'L') 1792 @:ASSERT(side == 'r' .or. side == 'R' .or. side == 'l' .or. side == 'L') 1793 1794 if (present(n)) then 1795 in = n 1796 else 1797 in = size(C,dim=2) 1798 end if 1799 if (present(m)) then 1800 im = m 1801 else 1802 im = size(C,dim=1) 1803 end if 1804 if (iUplo=='l'.or.iUplo=='L') then 1805 ka = im 1806 else 1807 ka = in 1808 end if 1809 1810 @:ASSERT(im>0) 1811 @:ASSERT(in>0) 1812 @:ASSERT((lda>=im).and.(iUplo=='l'.or.iUplo=='L').or.(lda>=in)) 1813 @:ASSERT(size(a,dim=2)>=ka) 1814 @:ASSERT(ldb>=im) 1815 @:ASSERT(ldc>=im) 1816 @:ASSERT(size(B,dim=2)>=in) 1817 @:ASSERT(size(C,dim=2)>=in) 1818 1819 call chemm(side, iUplo, im, in, iAlpha, A, lda, B, ldb, iBeta, C, ldc) 1820 1821 end subroutine hemm_cmplx 1822 1823 1824 !> double precision hermitian matrix * general matrix multiply 1825 subroutine hemm_dblecmplx(C, side, A, B, uplo, alpha, beta, m, n) 1826 1827 !> general matrix output 1828 complex(rdp), intent(inout) :: C(:,:) 1829 1830 !> symmetric matrix on 'l'eft or 'r'ight , where A is symmetric and B is gen 1831 !> 'l' C := alpha*A*B + beta*C, SIDE = 'R' or 'r' C := alpha*B*A + beta*C 1832 character, intent(in) :: side 1833 1834 !> hermitian matrix 1835 complex(rdp), intent(in) :: A(:,:) 1836 1837 !> general matrix 1838 complex(rdp), intent(in) :: B(:,:) 1839 1840 !> A is an 'U'pper or 'L'ower triangle matrix, defaults to lower 1841 character, intent(in), optional :: uplo 1842 1843 !> defaults to 1 if not set 1844 complex(rdp), intent(in), optional :: alpha 1845 1846 !> defaults to 0 if not set 1847 complex(rdp), intent(in), optional :: beta 1848 1849 !> specifies the number of rows of the matrix C 1850 integer, intent(in), optional :: m 1851 1852 !> specifies the number of columns of the matrix C 1853 integer, intent(in), optional :: n 1854 1855 integer :: lda, ldb, ldc, ka, im, in 1856 character :: iUplo 1857 complex(rdp) :: iAlpha, iBeta 1858 1859 if (present(uplo)) then 1860 iUplo = uplo 1861 else 1862 iUplo = 'L' 1863 end if 1864 if (present(alpha)) then 1865 iAlpha = alpha 1866 else 1867 iAlpha = (1.0_rdp, 0.0_rdp) 1868 end if 1869 if (present(beta)) then 1870 iBeta = beta 1871 else 1872 iBeta = (0.0_rdp, 0.0_rdp) 1873 end if 1874 1875 lda = size(a,dim=1) 1876 ldb = size(b,dim=1) 1877 ldc = size(c,dim=1) 1878 1879 @:ASSERT(iUplo == 'u' .or. iUplo == 'U' .or. iUplo == 'l' .or. iUplo == 'L') 1880 @:ASSERT(side == 'r' .or. side == 'R' .or. side == 'l' .or. side == 'L') 1881 1882 if (present(n)) then 1883 in = n 1884 else 1885 in = size(C,dim=2) 1886 end if 1887 if (present(m)) then 1888 im = m 1889 else 1890 im = size(C,dim=1) 1891 end if 1892 if (iUplo=='l'.or.iUplo=='L') then 1893 ka = im 1894 else 1895 ka = in 1896 end if 1897 1898 @:ASSERT(im>0) 1899 @:ASSERT(in>0) 1900 @:ASSERT((lda>=im).and.(iUplo=='l'.or.iUplo=='L').or.(lda>=in)) 1901 @:ASSERT(size(a,dim=2)>=ka) 1902 @:ASSERT(ldb>=im) 1903 @:ASSERT(ldc>=im) 1904 @:ASSERT(size(B,dim=2)>=in) 1905 @:ASSERT(size(C,dim=2)>=in) 1906 1907 call zhemm(side, iUplo, im, in, iAlpha, A, lda, B, ldb, iBeta, C, ldc) 1908 1909 end subroutine hemm_dblecmplx 1910 1911end module dftbp_blasroutines 1912