linalg 1.8.2
A linear algebra library that provides a user-friendly interface to several BLAS and LAPACK routines.
Loading...
Searching...
No Matches
linalg_basic.f90
1! linalg_basic.f90
2
3submodule(linalg) linalg_basic
4 use blas
5 use lapack
6 implicit none
7contains
8! ******************************************************************************
9! MATRIX MULTIPLICATION ROUTINES
10! ------------------------------------------------------------------------------
11 module subroutine mtx_mult_mtx(transa, transb, alpha, a, b, beta, c, err)
12 ! Arguments
13 logical, intent(in) :: transa, transb
14 real(real64), intent(in) :: alpha, beta
15 real(real64), intent(in), dimension(:,:) :: a, b
16 real(real64), intent(inout), dimension(:,:) :: c
17 class(errors), intent(inout), optional, target :: err
18
19 ! Parameters
20 real(real64), parameter :: zero = 0.0d0
21 real(real64), parameter :: one = 1.0d0
22
23 ! Local Variables
24 character :: ta, tb
25 integer(int32) :: m, n, k, lda, ldb, flag
26 class(errors), pointer :: errmgr
27 type(errors), target :: deferr
28 character(len = :), allocatable :: errmsg
29
30 ! Initialization
31 m = size(c, 1)
32 n = size(c, 2)
33 if (transa) then ! K = # of columns in op(A) (# of rows in op(B))
34 k = size(a, 1)
35 ta = 'T'
36 lda = k
37 else
38 k = size(a, 2)
39 ta = 'N'
40 lda = m
41 end if
42 if (transb) then
43 tb = 'T'
44 ldb = n
45 else
46 tb = 'N'
47 ldb = k
48 end if
49 if (present(err)) then
50 errmgr => err
51 else
52 errmgr => deferr
53 end if
54
55 ! Input Check
56 flag = 0
57 if (transa) then
58 if (size(a, 2) /= m) flag = 4
59 else
60 if (size(a, 1) /= m) flag = 4
61 end if
62 if (transb) then
63 if (size(b, 2) /= k .or. size(b, 1) /= n) flag = 5
64 else
65 if (size(b, 1) /= k .or. size(b, 2) /= n) flag = 5
66 end if
67 if (flag /= 0) then
68 ! ERROR: Matrix dimensions mismatch
69 allocate(character(len = 256) :: errmsg)
70 write(errmsg, 100) &
71 "Matrix dimension mismatch. Input number ", flag, &
72 " was not sized correctly."
73 call errmgr%report_error("mtx_mult_mtx", errmsg, &
74 la_array_size_error)
75 return
76 end if
77
78 ! Call DGEMM
79 call dgemm(ta, tb, m, n, k, alpha, a, lda, b, ldb, beta, c, m)
80
81 ! Formatting
82100 format(a, i0, a)
83 end subroutine
84
85! ------------------------------------------------------------------------------
86 module subroutine mtx_mult_vec(trans, alpha, a, b, beta, c, err)
87 ! Arguments
88 logical, intent(in) :: trans
89 real(real64), intent(in) :: alpha, beta
90 real(real64), intent(in), dimension(:,:) :: a
91 real(real64), intent(in), dimension(:) :: b
92 real(real64), intent(inout), dimension(:) :: c
93 class(errors), intent(inout), optional, target :: err
94
95 ! Local Variables
96 character :: t
97 integer(int32) :: m, n, flag
98 class(errors), pointer :: errmgr
99 type(errors), target :: deferr
100 character(len = :), allocatable :: errmsg
101
102 ! Initialization
103 m = size(a, 1)
104 n = size(a, 2)
105 t = 'N'
106 if (trans) t = 'T'
107 if (present(err)) then
108 errmgr => err
109 else
110 errmgr => deferr
111 end if
112
113 ! Input Check
114 flag = 0
115 if (trans) then
116 if (size(b) /= m) then
117 flag = 4
118 else if (size(c) /= n) then
119 flag = 6
120 end if
121 else
122 if (size(b) /= n) then
123 flag = 4
124 else if (size(c) /= m) then
125 flag = 6
126 end if
127 end if
128 if (flag /= 0) then
129 ! ERROR: Matrix dimensions mismatch
130 allocate(character(len = 256) :: errmsg)
131 write(errmsg, 100) &
132 "Matrix dimension mismatch. Input number ", flag, &
133 " was not sized correctly."
134 call errmgr%report_error("mtx_mult_vec", errmsg, &
135 la_array_size_error)
136 return
137 end if
138
139 ! Call DGEMV
140 call dgemv(t, m, n, alpha, a, m, b, 1, beta, c, 1)
141
142 ! Formatting
143100 format(a, i0, a)
144 end subroutine
145
146! xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx !
147! COMPLEX VALUED VERSIONS !
148! xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx !
149 module subroutine cmtx_mult_mtx(opa, opb, alpha, a, b, beta, c, err)
150 ! Arguments
151 integer(int32), intent(in) :: opa, opb
152 complex(real64), intent(in) :: alpha, beta
153 complex(real64), intent(in), dimension(:,:) :: a, b
154 complex(real64), intent(inout), dimension(:,:) :: c
155 class(errors), intent(inout), optional, target :: err
156
157 ! Parameters
158 real(real64), parameter :: zero = 0.0d0
159 real(real64), parameter :: one = 1.0d0
160
161 ! Local Variables
162 character :: ta, tb
163 integer(int32) :: m, n, k, lda, ldb, flag
164 class(errors), pointer :: errmgr
165 type(errors), target :: deferr
166 character(len = :), allocatable :: errmsg
167
168 ! Initialization
169 m = size(c, 1)
170 n = size(c, 2)
171 if (opa == la_transpose) then ! K = # of columns in op(A) (# of rows in op(B))
172 k = size(a, 1)
173 ta = 'T'
174 lda = k
175 else if (opa == la_hermitian_transpose) then
176 k = size(a, 1)
177 ta = 'C'
178 lda = k
179 else
180 k = size(a, 2)
181 ta = 'N'
182 lda = m
183 end if
184 if (opb == la_transpose) then
185 tb = 'T'
186 ldb = n
187 else if (opb == la_hermitian_transpose) then
188 tb = 'C'
189 ldb = n
190 else
191 tb = 'N'
192 ldb = k
193 end if
194 if (present(err)) then
195 errmgr => err
196 else
197 errmgr => deferr
198 end if
199
200 ! Input Check
201 flag = 0
202 if (opa == la_transpose .or. opa == la_hermitian_transpose) then
203 if (size(a, 2) /= m) flag = 4
204 else
205 if (size(a, 1) /= m) flag = 4
206 end if
207 if (opb == la_transpose .or. opb == la_hermitian_transpose) then
208 if (size(b, 2) /= k .or. size(b, 1) /= n) flag = 5
209 else
210 if (size(b, 1) /= k .or. size(b, 2) /= n) flag = 5
211 end if
212 if (flag /= 0) then
213 ! ERROR: Matrix dimensions mismatch
214 allocate(character(len = 256) :: errmsg)
215 write(errmsg, 100) &
216 "Matrix dimension mismatch. Input number ", flag, &
217 " was not sized correctly."
218 call errmgr%report_error("cmtx_mult_mtx", errmsg, &
219 la_array_size_error)
220 return
221 end if
222
223 ! Call ZGEMM
224 call zgemm(ta, tb, m, n, k, alpha, a, lda, b, ldb, beta, c, m)
225
226 ! Formatting
227100 format(a, i0, a)
228 end subroutine
229
230! ------------------------------------------------------------------------------
231 module subroutine cmtx_mult_vec(opa, alpha, a, b, beta, c, err)
232 ! Arguments
233 integer(int32), intent(in) :: opa
234 complex(real64), intent(in) :: alpha, beta
235 complex(real64), intent(in), dimension(:,:) :: a
236 complex(real64), intent(in), dimension(:) :: b
237 complex(real64), intent(inout), dimension(:) :: c
238 class(errors), intent(inout), optional, target :: err
239
240 ! Local Variables
241 character :: t
242 integer(int32) :: m, n, flag
243 class(errors), pointer :: errmgr
244 type(errors), target :: deferr
245 character(len = :), allocatable :: errmsg
246
247 ! Initialization
248 m = size(a, 1)
249 n = size(a, 2)
250 if (opa == la_transpose) then
251 t = 'T'
252 else if (opa == la_hermitian_transpose) then
253 t = 'C'
254 else
255 t = 'N'
256 end if
257 if (present(err)) then
258 errmgr => err
259 else
260 errmgr => deferr
261 end if
262
263 ! Input Check
264 flag = 0
265 if (opa == la_transpose .or. opa == la_hermitian_transpose) then
266 if (size(b) /= m) then
267 flag = 4
268 else if (size(c) /= n) then
269 flag = 6
270 end if
271 else
272 if (size(b) /= n) then
273 flag = 4
274 else if (size(c) /= m) then
275 flag = 6
276 end if
277 end if
278 if (flag /= 0) then
279 ! ERROR: Matrix dimensions mismatch
280 allocate(character(len = 256) :: errmsg)
281 write(errmsg, 100) &
282 "Matrix dimension mismatch. Input number ", flag, &
283 " was not sized correctly."
284 call errmgr%report_error("cmtx_mult_vec", errmsg, &
285 la_array_size_error)
286 return
287 end if
288
289 ! Call ZGEMV
290 call zgemv(t, m, n, alpha, a, m, b, 1, beta, c, 1)
291
292 ! Formatting
293100 format(a, i0, a)
294 end subroutine
295
296! ******************************************************************************
297! RANK 1 UPDATE
298! ------------------------------------------------------------------------------
299 module subroutine rank1_update_dbl(alpha, x, y, a, err)
300 ! Arguments
301 real(real64), intent(in) :: alpha
302 real(real64), intent(in), dimension(:) :: x, y
303 real(real64), intent(inout), dimension(:,:) :: a
304 class(errors), intent(inout), optional, target :: err
305
306 ! Parameters
307 real(real64), parameter :: zero = 0.0d0
308
309 ! Local Variables
310 integer(int32) :: j, m, n
311 real(real64) :: temp
312 class(errors), pointer :: errmgr
313 type(errors), target :: deferr
314
315 ! Initialization
316 m = size(x)
317 n = size(y)
318 if (present(err)) then
319 errmgr => err
320 else
321 errmgr => deferr
322 end if
323
324 ! Input Check
325 if (size(a, 1) /= m .or. size(a, 2) /= n) then
326 ! ERROR: Matrix dimension array
327 call errmgr%report_error("rank1_update_dbl", &
328 "Matrix dimension mismatch.", la_array_size_error)
329 return
330 end if
331
332 ! Process
333 do j = 1, n
334 if (y(j) /= zero) then
335 temp = alpha * y(j)
336 a(:,j) = a(:,j) + temp * x
337 end if
338 end do
339 end subroutine
340
341! xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx !
342! COMPLEX VALUED VERSIONS !
343! xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx !
344 module subroutine rank1_update_cmplx(alpha, x, y, a, err)
345 ! Arguments
346 complex(real64), intent(in) :: alpha
347 complex(real64), intent(in), dimension(:) :: x, y
348 complex(real64), intent(inout), dimension(:,:) :: a
349 class(errors), intent(inout), optional, target :: err
350
351 ! Parameters
352 complex(real64), parameter :: zero = (0.0d0, 0.0d0)
353
354 ! Local Variables
355 integer(int32) :: j, m, n
356 complex(real64) :: temp
357 class(errors), pointer :: errmgr
358 type(errors), target :: deferr
359
360 ! Initialization
361 m = size(x)
362 n = size(y)
363 if (present(err)) then
364 errmgr => err
365 else
366 errmgr => deferr
367 end if
368
369 ! Input Check
370 if (size(a, 1) /= m .or. size(a, 2) /= n) then
371 ! ERROR: Matrix dimension array
372 call errmgr%report_error("rank1_update_cmplx", &
373 "Matrix dimension mismatch.", la_array_size_error)
374 return
375 end if
376
377 ! Process
378 do j = 1, n
379 if (y(j) /= zero) then
380 temp = alpha * conjg(y(j))
381 a(:,j) = a(:,j) + temp * x
382 end if
383 end do
384 end subroutine
385
386! ******************************************************************************
387! DIAGONAL MATRIX MULTIPLICATION ROUTINES
388! ------------------------------------------------------------------------------
389 module subroutine diag_mtx_mult_mtx(lside, trans, alpha, a, b, beta, c, err)
390 ! Arguments
391 logical, intent(in) :: lside, trans
392 real(real64) :: alpha, beta
393 real(real64), intent(in), dimension(:) :: a
394 real(real64), intent(in), dimension(:,:) :: b
395 real(real64), intent(inout), dimension(:,:) :: c
396 class(errors), intent(inout), optional, target :: err
397
398 ! Parameters
399 real(real64), parameter :: zero = 0.0d0
400 real(real64), parameter :: one = 1.0d0
401
402 ! Local Variables
403 integer(int32) :: i, m, n, k, nrowb, ncolb, flag
404 real(real64) :: temp
405 class(errors), pointer :: errmgr
406 type(errors), target :: deferr
407 character(len = :), allocatable :: errmsg
408
409 ! Initialization
410 m = size(c, 1)
411 n = size(c, 2)
412 k = size(a)
413 nrowb = size(b, 1)
414 ncolb = size(b, 2)
415 if (present(err)) then
416 errmgr => err
417 else
418 errmgr => deferr
419 end if
420
421 ! Input Check
422 flag = 0
423 if (lside) then
424 if (k > m) then
425 flag = 4
426 else
427 if (trans) then
428 ! Compute C = alpha * A * B**T + beta * C
429 if (nrowb /= n .or. ncolb < k) flag = 5
430 else
431 ! Compute C = alpha * A * B + beta * C
432 if (nrowb < k .or. ncolb /= n) flag = 5
433 end if
434 end if
435 else
436 if (k > n) then
437 flag = 4
438 else
439 if (trans) then
440 ! Compute C = alpha * B**T * A + beta * C
441 if (ncolb /= m .or. nrowb < k) flag = 5
442 else
443 ! Compute C = alpha * B * A + beta * C
444 if (nrowb /= m .or. ncolb < k) flag = 5
445 end if
446 end if
447 end if
448 if (flag /= 0) then
449 ! ERROR: One of the input arrays is not sized correctly
450 allocate(character(len = 256) :: errmsg)
451 write(errmsg, 100) "Input number ", flag, &
452 " is not sized correctly."
453 call errmgr%report_error("diag_mtx_mult_mtx", trim(errmsg), &
454 la_array_size_error)
455 return
456 end if
457
458 ! Deal with ALPHA == 0
459 if (alpha == 0) then
460 if (beta == zero) then
461 c = zero
462 else if (beta /= one) then
463 c = beta * c
464 end if
465 return
466 end if
467
468 ! Process
469 if (lside) then
470 if (trans) then
471 ! Compute C = alpha * A * B**T + beta * C
472 do i = 1, k
473 if (beta == zero) then
474 c(i,:) = zero
475 else if (beta /= one) then
476 c(i,:) = beta * c(i,:)
477 end if
478 temp = alpha * a(i)
479 c(i,:) = c(i,:) + temp * b(:,i)
480 end do
481 else
482 ! Compute C = alpha * A * B + beta * C
483 do i = 1, k
484 if (beta == zero) then
485 c(i,:) = zero
486 else if (beta /= one) then
487 c(i,:) = beta * c(i,:)
488 end if
489 temp = alpha * a(i)
490 c(i,:) = c(i,:) + temp * b(i,:)
491 end do
492 end if
493
494 ! Handle extra rows
495 if (m > k) then
496 if (beta == zero) then
497 c(k+1:m,:) = zero
498 else
499 c(k+1:m,:) = beta * c(k+1:m,:)
500 end if
501 end if
502 else
503 if (trans) then
504 ! Compute C = alpha * B**T * A + beta * C
505 do i = 1, k
506 if (beta == zero) then
507 c(:,i) = zero
508 else if (beta /= one) then
509 c(:,i) = beta * c(:,i)
510 end if
511 temp = alpha * a(i)
512 c(:,i) = c(:,i) + temp * b(i,:)
513 end do
514 else
515 ! Compute C = alpha * B * A + beta * C
516 do i = 1, k
517 if (beta == zero) then
518 c(:,i) = zero
519 else if (beta /= one) then
520 c(:,i) = beta * c(:,i)
521 end if
522 temp = alpha * a(i)
523 c(:,i) = c(:,i) + temp * b(:,i)
524 end do
525 end if
526
527 ! Handle extra columns
528 if (n > k) then
529 if (beta == zero) then
530 c(:,k+1:m) = zero
531 else if (beta /= one) then
532 c(:,k+1:m) = beta * c(:,k+1:m)
533 end if
534 end if
535 end if
536
537 ! Formatting
538100 format(a, i0, a)
539 end subroutine
540
541! ------------------------------------------------------------------------------
542 module subroutine diag_mtx_mult_mtx2(lside, alpha, a, b, err)
543 ! Arguments
544 logical, intent(in) :: lside
545 real(real64), intent(in) :: alpha
546 real(real64), intent(in), dimension(:) :: a
547 real(real64), intent(inout), dimension(:,:) :: b
548 class(errors), intent(inout), optional, target :: err
549
550 ! Parameters
551 real(real64), parameter :: zero = 0.0d0
552 real(real64), parameter :: one = 1.0d0
553
554 ! Local Variables
555 integer(int32) :: i, m, n, k
556 real(real64) :: temp
557 class(errors), pointer :: errmgr
558 type(errors), target :: deferr
559
560 ! Initialization
561 m = size(b, 1)
562 n = size(b, 2)
563 k = size(a)
564 if (present(err)) then
565 errmgr => err
566 else
567 errmgr => deferr
568 end if
569
570 ! Input Check
571 if ((lside .and. k > m) .or. (.not.lside .and. k > n)) then
572 ! ERROR: One of the input arrays is not sized correctly
573 call errmgr%report_error("diag_mtx_mult_mtx2", &
574 "Input number 3 is not sized correctly.", &
575 la_array_size_error)
576 return
577 end if
578
579 ! Process
580 if (lside) then
581 ! Compute B = alpha * A * B
582 do i = 1, k
583 temp = alpha * a(i)
584 b(i,:) = temp * b(i,:)
585 end do
586 if (m > k) b(k+1:m,:) = zero
587 else
588 ! Compute B = alpha * B * A
589 do i = 1, k
590 temp = alpha * a(i)
591 b(:,i) = temp * b(:,i)
592 end do
593 if (n > k) b(:,k+1:n) = zero
594 end if
595 end subroutine
596
597! ------------------------------------------------------------------------------
598 module subroutine diag_mtx_mult_mtx3(lside, trans, alpha, a, b, beta, c, err)
599 ! Arguments
600 logical, intent(in) :: lside, trans
601 real(real64) :: alpha, beta
602 complex(real64), intent(in), dimension(:) :: a
603 real(real64), intent(in), dimension(:,:) :: b
604 complex(real64), intent(inout), dimension(:,:) :: c
605 class(errors), intent(inout), optional, target :: err
606
607 ! Parameters
608 complex(real64), parameter :: zero = (0.0d0, 0.0d0)
609 complex(real64), parameter :: one = (1.0d0, 0.0d0)
610
611 ! Local Variables
612 integer(int32) :: i, m, n, k, nrowb, ncolb, flag
613 complex(real64) :: temp
614 class(errors), pointer :: errmgr
615 type(errors), target :: deferr
616 character(len = :), allocatable :: errmsg
617
618 ! Initialization
619 m = size(c, 1)
620 n = size(c, 2)
621 k = size(a)
622 nrowb = size(b, 1)
623 ncolb = size(b, 2)
624 if (present(err)) then
625 errmgr => err
626 else
627 errmgr => deferr
628 end if
629
630 ! Input Check
631 flag = 0
632 if (lside) then
633 if (k > m) then
634 flag = 4
635 else
636 if (trans) then
637 ! Compute C = alpha * A * B**T + beta * C
638 if (nrowb /= n .or. ncolb < k) flag = 5
639 else
640 ! Compute C = alpha * A * B + beta * C
641 if (nrowb < k .or. ncolb /= n) flag = 5
642 end if
643 end if
644 else
645 if (k > n) then
646 flag = 4
647 else
648 if (trans) then
649 ! Compute C = alpha * B**T * A + beta * C
650 if (ncolb /= m .or. nrowb < k) flag = 5
651 else
652 ! Compute C = alpha * B * A + beta * C
653 if (nrowb /= m .or. ncolb < k) flag = 5
654 end if
655 end if
656 end if
657 if (flag /= 0) then
658 ! ERROR: One of the input arrays is not sized correctly
659 allocate(character(len = 256) :: errmsg)
660 write(errmsg, 100) "Input number ", flag, &
661 " is not sized correctly."
662 call errmgr%report_error("diag_mtx_mult_mtx3", trim(errmsg), &
663 la_array_size_error)
664 return
665 end if
666
667 ! Deal with ALPHA == 0
668 if (alpha == 0) then
669 if (beta == zero) then
670 c = zero
671 else if (beta /= one) then
672 c = beta * c
673 end if
674 return
675 end if
676
677 ! Process
678 if (lside) then
679 if (trans) then
680 ! Compute C = alpha * A * B**T + beta * C
681 do i = 1, k
682 if (beta == zero) then
683 c(i,:) = zero
684 else if (beta /= one) then
685 c(i,:) = beta * c(i,:)
686 end if
687 temp = alpha * a(i)
688 c(i,:) = c(i,:) + temp * b(:,i)
689 end do
690 else
691 ! Compute C = alpha * A * B + beta * C
692 do i = 1, k
693 if (beta == zero) then
694 c(i,:) = zero
695 else if (beta /= one) then
696 c(i,:) = beta * c(i,:)
697 end if
698 temp = alpha * a(i)
699 c(i,:) = c(i,:) + temp * b(i,:)
700 end do
701 end if
702
703 ! Handle extra rows
704 if (m > k) then
705 if (beta == zero) then
706 c(k+1:m,:) = zero
707 else
708 c(k+1:m,:) = beta * c(k+1:m,:)
709 end if
710 end if
711 else
712 if (trans) then
713 ! Compute C = alpha * B**T * A + beta * C
714 do i = 1, k
715 if (beta == zero) then
716 c(:,i) = zero
717 else if (beta /= one) then
718 c(:,i) = beta * c(:,i)
719 end if
720 temp = alpha * a(i)
721 c(:,i) = c(:,i) + temp * b(i,:)
722 end do
723 else
724 ! Compute C = alpha * B * A + beta * C
725 do i = 1, k
726 if (beta == zero) then
727 c(:,i) = zero
728 else if (beta /= one) then
729 c(:,i) = beta * c(:,i)
730 end if
731 temp = alpha * a(i)
732 c(:,i) = c(:,i) + temp * b(:,i)
733 end do
734 end if
735
736 ! Handle extra columns
737 if (n > k) then
738 if (beta == zero) then
739 c(:,k+1:m) = zero
740 else if (beta /= one) then
741 c(:,k+1:m) = beta * c(:,k+1:m)
742 end if
743 end if
744 end if
745
746 ! Formatting
747100 format(a, i0, a)
748 end subroutine
749
750! ------------------------------------------------------------------------------
751 module subroutine diag_mtx_mult_mtx4(lside, opb, alpha, a, b, beta, c, err)
752 ! Arguments
753 logical, intent(in) :: lside
754 integer(int32), intent(in) :: opb
755 real(real64) :: alpha, beta
756 complex(real64), intent(in), dimension(:) :: a
757 complex(real64), intent(in), dimension(:,:) :: b
758 complex(real64), intent(inout), dimension(:,:) :: c
759 class(errors), intent(inout), optional, target :: err
760
761 ! Parameters
762 complex(real64), parameter :: zero = (0.0d0, 0.0d0)
763 complex(real64), parameter :: one = (1.0d0, 0.0d0)
764
765 ! Local Variables
766 integer(int32) :: i, m, n, k, nrowb, ncolb, flag
767 complex(real64) :: temp
768 class(errors), pointer :: errmgr
769 type(errors), target :: deferr
770 character(len = :), allocatable :: errmsg
771
772 ! Initialization
773 m = size(c, 1)
774 n = size(c, 2)
775 k = size(a)
776 nrowb = size(b, 1)
777 ncolb = size(b, 2)
778 if (present(err)) then
779 errmgr => err
780 else
781 errmgr => deferr
782 end if
783
784 ! Input Check
785 flag = 0
786 if (lside) then
787 if (k > m) then
788 flag = 4
789 else
790 if (opb == la_transpose .or. opb == la_hermitian_transpose) then
791 ! Compute C = alpha * A * B**T + beta * C
792 if (nrowb /= n .or. ncolb < k) flag = 5
793 else
794 ! Compute C = alpha * A * B + beta * C
795 if (nrowb < k .or. ncolb /= n) flag = 5
796 end if
797 end if
798 else
799 if (k > n) then
800 flag = 4
801 else
802 if (opb == la_transpose .or. opb == la_hermitian_transpose) then
803 ! Compute C = alpha * B**T * A + beta * C
804 if (ncolb /= m .or. nrowb < k) flag = 5
805 else
806 ! Compute C = alpha * B * A + beta * C
807 if (nrowb /= m .or. ncolb < k) flag = 5
808 end if
809 end if
810 end if
811 if (flag /= 0) then
812 ! ERROR: One of the input arrays is not sized correctly
813 allocate(character(len = 256) :: errmsg)
814 write(errmsg, 100) "Input number ", flag, &
815 " is not sized correctly."
816 call errmgr%report_error("diag_mtx_mult_mtx4", trim(errmsg), &
817 la_array_size_error)
818 return
819 end if
820
821 ! Deal with ALPHA == 0
822 if (alpha == 0) then
823 if (beta == zero) then
824 c = zero
825 else if (beta /= one) then
826 c = beta * c
827 end if
828 return
829 end if
830
831 ! Process
832 if (lside) then
833 if (opb == la_transpose) then
834 ! Compute C = alpha * A * B**T + beta * C
835 do i = 1, k
836 if (beta == zero) then
837 c(i,:) = zero
838 else if (beta /= one) then
839 c(i,:) = beta * c(i,:)
840 end if
841 temp = alpha * a(i)
842 c(i,:) = c(i,:) + temp * b(:,i)
843 end do
844 else if (opb == la_hermitian_transpose) then
845 ! Compute C = alpha * A * B**H + beta * C
846 do i = 1, k
847 if (beta == zero) then
848 c(i,:) = zero
849 else if (beta /= one) then
850 c(i,:) = beta * c(i,:)
851 end if
852 temp = alpha * a(i)
853 c(i,:) = c(i,:) + temp * conjg(b(:,i))
854 end do
855 else
856 ! Compute C = alpha * A * B + beta * C
857 do i = 1, k
858 if (beta == zero) then
859 c(i,:) = zero
860 else if (beta /= one) then
861 c(i,:) = beta * c(i,:)
862 end if
863 temp = alpha * a(i)
864 c(i,:) = c(i,:) + temp * b(i,:)
865 end do
866 end if
867
868 ! Handle extra rows
869 if (m > k) then
870 if (beta == zero) then
871 c(k+1:m,:) = zero
872 else
873 c(k+1:m,:) = beta * c(k+1:m,:)
874 end if
875 end if
876 else
877 if (opb == la_transpose) then
878 ! Compute C = alpha * B**T * A + beta * C
879 do i = 1, k
880 if (beta == zero) then
881 c(:,i) = zero
882 else if (beta /= one) then
883 c(:,i) = beta * c(:,i)
884 end if
885 temp = alpha * a(i)
886 c(:,i) = c(:,i) + temp * b(i,:)
887 end do
888 else if (opb == la_hermitian_transpose) then
889 ! Compute C = alpha * B**H * A + beta * C
890 do i = 1, k
891 if (beta == zero) then
892 c(:,i) = zero
893 else if (beta /= one) then
894 c(:,i) = beta * c(:,i)
895 end if
896 temp = alpha * a(i)
897 c(:,i) = c(:,i) + temp * conjg(b(i,:))
898 end do
899 else
900 ! Compute C = alpha * B * A + beta * C
901 do i = 1, k
902 if (beta == zero) then
903 c(:,i) = zero
904 else if (beta /= one) then
905 c(:,i) = beta * c(:,i)
906 end if
907 temp = alpha * a(i)
908 c(:,i) = c(:,i) + temp * b(:,i)
909 end do
910 end if
911
912 ! Handle extra columns
913 if (n > k) then
914 if (beta == zero) then
915 c(:,k+1:m) = zero
916 else if (beta /= one) then
917 c(:,k+1:m) = beta * c(:,k+1:m)
918 end if
919 end if
920 end if
921
922 ! Formatting
923100 format(a, i0, a)
924 end subroutine
925
926! ------------------------------------------------------------------------------
927 module subroutine diag_mtx_mult_mtx_cmplx(lside, opb, alpha, a, b, beta, c, err)
928 ! Arguments
929 logical, intent(in) :: lside
930 integer(int32), intent(in) :: opb
931 complex(real64) :: alpha, beta
932 complex(real64), intent(in), dimension(:) :: a
933 complex(real64), intent(in), dimension(:,:) :: b
934 complex(real64), intent(inout), dimension(:,:) :: c
935 class(errors), intent(inout), optional, target :: err
936
937 ! Parameters
938 complex(real64), parameter :: zero = (0.0d0, 0.0d0)
939 complex(real64), parameter :: one = (1.0d0, 0.0d0)
940
941 ! Local Variables
942 integer(int32) :: i, m, n, k, nrowb, ncolb, flag
943 complex(real64) :: temp
944 class(errors), pointer :: errmgr
945 type(errors), target :: deferr
946 character(len = :), allocatable :: errmsg
947
948 ! Initialization
949 m = size(c, 1)
950 n = size(c, 2)
951 k = size(a)
952 nrowb = size(b, 1)
953 ncolb = size(b, 2)
954 if (present(err)) then
955 errmgr => err
956 else
957 errmgr => deferr
958 end if
959
960 ! Input Check
961 flag = 0
962 if (lside) then
963 if (k > m) then
964 flag = 4
965 else
966 if (opb == la_transpose .or. opb == la_hermitian_transpose) then
967 ! Compute C = alpha * A * B**T + beta * C
968 if (nrowb /= n .or. ncolb < k) flag = 5
969 else
970 ! Compute C = alpha * A * B + beta * C
971 if (nrowb < k .or. ncolb /= n) flag = 5
972 end if
973 end if
974 else
975 if (k > n) then
976 flag = 4
977 else
978 if (opb == la_transpose .or. opb == la_hermitian_transpose) then
979 ! Compute C = alpha * B**T * A + beta * C
980 if (ncolb /= m .or. nrowb < k) flag = 5
981 else
982 ! Compute C = alpha * B * A + beta * C
983 if (nrowb /= m .or. ncolb < k) flag = 5
984 end if
985 end if
986 end if
987 if (flag /= 0) then
988 ! ERROR: One of the input arrays is not sized correctly
989 allocate(character(len = 256) :: errmsg)
990 write(errmsg, 100) "Input number ", flag, &
991 " is not sized correctly."
992 call errmgr%report_error("diag_mtx_mult_mtx_cmplx", trim(errmsg), &
993 la_array_size_error)
994 return
995 end if
996
997 ! Deal with ALPHA == 0
998 if (alpha == 0) then
999 if (beta == zero) then
1000 c = zero
1001 else if (beta /= one) then
1002 c = beta * c
1003 end if
1004 return
1005 end if
1006
1007 ! Process
1008 if (lside) then
1009 if (opb == la_transpose) then
1010 ! Compute C = alpha * A * B**T + beta * C
1011 do i = 1, k
1012 if (beta == zero) then
1013 c(i,:) = zero
1014 else if (beta /= one) then
1015 c(i,:) = beta * c(i,:)
1016 end if
1017 temp = alpha * a(i)
1018 c(i,:) = c(i,:) + temp * b(:,i)
1019 end do
1020 else if (opb == la_hermitian_transpose) then
1021 ! Compute C = alpha * A * B**H + beta * C
1022 do i = 1, k
1023 if (beta == zero) then
1024 c(i,:) = zero
1025 else if (beta /= one) then
1026 c(i,:) = beta * c(i,:)
1027 end if
1028 temp = alpha * a(i)
1029 c(i,:) = c(i,:) + temp * conjg(b(:,i))
1030 end do
1031 else
1032 ! Compute C = alpha * A * B + beta * C
1033 do i = 1, k
1034 if (beta == zero) then
1035 c(i,:) = zero
1036 else if (beta /= one) then
1037 c(i,:) = beta * c(i,:)
1038 end if
1039 temp = alpha * a(i)
1040 c(i,:) = c(i,:) + temp * b(i,:)
1041 end do
1042 end if
1043
1044 ! Handle extra rows
1045 if (m > k) then
1046 if (beta == zero) then
1047 c(k+1:m,:) = zero
1048 else
1049 c(k+1:m,:) = beta * c(k+1:m,:)
1050 end if
1051 end if
1052 else
1053 if (opb == la_transpose) then
1054 ! Compute C = alpha * B**T * A + beta * C
1055 do i = 1, k
1056 if (beta == zero) then
1057 c(:,i) = zero
1058 else if (beta /= one) then
1059 c(:,i) = beta * c(:,i)
1060 end if
1061 temp = alpha * a(i)
1062 c(:,i) = c(:,i) + temp * b(i,:)
1063 end do
1064 else if (opb == la_hermitian_transpose) then
1065 ! Compute C = alpha * B**H * A + beta * C
1066 do i = 1, k
1067 if (beta == zero) then
1068 c(:,i) = zero
1069 else if (beta /= one) then
1070 c(:,i) = beta * c(:,i)
1071 end if
1072 temp = alpha * a(i)
1073 c(:,i) = c(:,i) + temp * conjg(b(i,:))
1074 end do
1075 else
1076 ! Compute C = alpha * B * A + beta * C
1077 do i = 1, k
1078 if (beta == zero) then
1079 c(:,i) = zero
1080 else if (beta /= one) then
1081 c(:,i) = beta * c(:,i)
1082 end if
1083 temp = alpha * a(i)
1084 c(:,i) = c(:,i) + temp * b(:,i)
1085 end do
1086 end if
1087
1088 ! Handle extra columns
1089 if (n > k) then
1090 if (beta == zero) then
1091 c(:,k+1:m) = zero
1092 else if (beta /= one) then
1093 c(:,k+1:m) = beta * c(:,k+1:m)
1094 end if
1095 end if
1096 end if
1097
1098 ! Formatting
1099100 format(a, i0, a)
1100 end subroutine
1101
1102! ------------------------------------------------------------------------------
1103 module subroutine diag_mtx_mult_mtx2_cmplx(lside, alpha, a, b, err)
1104 ! Arguments
1105 logical, intent(in) :: lside
1106 complex(real64), intent(in) :: alpha
1107 complex(real64), intent(in), dimension(:) :: a
1108 complex(real64), intent(inout), dimension(:,:) :: b
1109 class(errors), intent(inout), optional, target :: err
1110
1111 ! Parameters
1112 complex(real64), parameter :: zero = (0.0d0, 0.0d0)
1113 complex(real64), parameter :: one = (1.0d0, 0.0d0)
1114
1115 ! Local Variables
1116 integer(int32) :: i, m, n, k
1117 complex(real64) :: temp
1118 class(errors), pointer :: errmgr
1119 type(errors), target :: deferr
1120
1121 ! Initialization
1122 m = size(b, 1)
1123 n = size(b, 2)
1124 k = size(a)
1125 if (present(err)) then
1126 errmgr => err
1127 else
1128 errmgr => deferr
1129 end if
1130
1131 ! Input Check
1132 if ((lside .and. k > m) .or. (.not.lside .and. k > n)) then
1133 ! ERROR: One of the input arrays is not sized correctly
1134 call errmgr%report_error("diag_mtx_mult_mtx2_cmplx", &
1135 "Input number 3 is not sized correctly.", &
1136 la_array_size_error)
1137 return
1138 end if
1139
1140 ! Process
1141 if (lside) then
1142 ! Compute B = alpha * A * B
1143 do i = 1, k
1144 temp = alpha * a(i)
1145 b(i,:) = temp * b(i,:)
1146 end do
1147 if (m > k) b(k+1:m,:) = zero
1148 else
1149 ! Compute B = alpha * B * A
1150 do i = 1, k
1151 temp = alpha * a(i)
1152 b(:,i) = temp * b(:,i)
1153 end do
1154 if (n > k) b(:,k+1:n) = zero
1155 end if
1156 end subroutine
1157
1158! ------------------------------------------------------------------------------
1159 module subroutine diag_mtx_mult_mtx_mix(lside, opb, alpha, a, b, beta, c, err)
1160 ! Arguments
1161 logical, intent(in) :: lside
1162 integer(int32), intent(in) :: opb
1163 complex(real64) :: alpha, beta
1164 real(real64), intent(in), dimension(:) :: a
1165 complex(real64), intent(in), dimension(:,:) :: b
1166 complex(real64), intent(inout), dimension(:,:) :: c
1167 class(errors), intent(inout), optional, target :: err
1168
1169 ! Parameters
1170 complex(real64), parameter :: zero = (0.0d0, 0.0d0)
1171 complex(real64), parameter :: one = (1.0d0, 0.0d0)
1172
1173 ! Local Variables
1174 integer(int32) :: i, m, n, k, nrowb, ncolb, flag
1175 complex(real64) :: temp
1176 class(errors), pointer :: errmgr
1177 type(errors), target :: deferr
1178 character(len = :), allocatable :: errmsg
1179
1180 ! Initialization
1181 m = size(c, 1)
1182 n = size(c, 2)
1183 k = size(a)
1184 nrowb = size(b, 1)
1185 ncolb = size(b, 2)
1186 if (present(err)) then
1187 errmgr => err
1188 else
1189 errmgr => deferr
1190 end if
1191
1192 ! Input Check
1193 flag = 0
1194 if (lside) then
1195 if (k > m) then
1196 flag = 4
1197 else
1198 if (opb == la_transpose .or. opb == la_hermitian_transpose) then
1199 ! Compute C = alpha * A * B**T + beta * C
1200 if (nrowb /= n .or. ncolb < k) flag = 5
1201 else
1202 ! Compute C = alpha * A * B + beta * C
1203 if (nrowb < k .or. ncolb /= n) flag = 5
1204 end if
1205 end if
1206 else
1207 if (k > n) then
1208 flag = 4
1209 else
1210 if (opb == la_transpose .or. opb == la_hermitian_transpose) then
1211 ! Compute C = alpha * B**T * A + beta * C
1212 if (ncolb /= m .or. nrowb < k) flag = 5
1213 else
1214 ! Compute C = alpha * B * A + beta * C
1215 if (nrowb /= m .or. ncolb < k) flag = 5
1216 end if
1217 end if
1218 end if
1219 if (flag /= 0) then
1220 ! ERROR: One of the input arrays is not sized correctly
1221 allocate(character(len = 256) :: errmsg)
1222 write(errmsg, 100) "Input number ", flag, &
1223 " is not sized correctly."
1224 call errmgr%report_error("diag_mtx_mult_mtx_mix", trim(errmsg), &
1225 la_array_size_error)
1226 return
1227 end if
1228
1229 ! Deal with ALPHA == 0
1230 if (alpha == 0) then
1231 if (beta == zero) then
1232 c = zero
1233 else if (beta /= one) then
1234 c = beta * c
1235 end if
1236 return
1237 end if
1238
1239 ! Process
1240 if (lside) then
1241 if (opb == la_transpose) then
1242 ! Compute C = alpha * A * B**T + beta * C
1243 do i = 1, k
1244 if (beta == zero) then
1245 c(i,:) = zero
1246 else if (beta /= one) then
1247 c(i,:) = beta * c(i,:)
1248 end if
1249 temp = alpha * a(i)
1250 c(i,:) = c(i,:) + temp * b(:,i)
1251 end do
1252 else if (opb == la_hermitian_transpose) then
1253 ! Compute C = alpha * A * B**H + beta * C
1254 do i = 1, k
1255 if (beta == zero) then
1256 c(i,:) = zero
1257 else if (beta /= one) then
1258 c(i,:) = beta * c(i,:)
1259 end if
1260 temp = alpha * a(i)
1261 c(i,:) = c(i,:) + temp * conjg(b(:,i))
1262 end do
1263 else
1264 ! Compute C = alpha * A * B + beta * C
1265 do i = 1, k
1266 if (beta == zero) then
1267 c(i,:) = zero
1268 else if (beta /= one) then
1269 c(i,:) = beta * c(i,:)
1270 end if
1271 temp = alpha * a(i)
1272 c(i,:) = c(i,:) + temp * b(i,:)
1273 end do
1274 end if
1275
1276 ! Handle extra rows
1277 if (m > k) then
1278 if (beta == zero) then
1279 c(k+1:m,:) = zero
1280 else
1281 c(k+1:m,:) = beta * c(k+1:m,:)
1282 end if
1283 end if
1284 else
1285 if (opb == la_transpose) then
1286 ! Compute C = alpha * B**T * A + beta * C
1287 do i = 1, k
1288 if (beta == zero) then
1289 c(:,i) = zero
1290 else if (beta /= one) then
1291 c(:,i) = beta * c(:,i)
1292 end if
1293 temp = alpha * a(i)
1294 c(:,i) = c(:,i) + temp * b(i,:)
1295 end do
1296 else if (opb == la_hermitian_transpose) then
1297 ! Compute C = alpha * B**H * A + beta * C
1298 do i = 1, k
1299 if (beta == zero) then
1300 c(:,i) = zero
1301 else if (beta /= one) then
1302 c(:,i) = beta * c(:,i)
1303 end if
1304 temp = alpha * a(i)
1305 c(:,i) = c(:,i) + temp * conjg(b(i,:))
1306 end do
1307 else
1308 ! Compute C = alpha * B * A + beta * C
1309 do i = 1, k
1310 if (beta == zero) then
1311 c(:,i) = zero
1312 else if (beta /= one) then
1313 c(:,i) = beta * c(:,i)
1314 end if
1315 temp = alpha * a(i)
1316 c(:,i) = c(:,i) + temp * b(:,i)
1317 end do
1318 end if
1319
1320 ! Handle extra columns
1321 if (n > k) then
1322 if (beta == zero) then
1323 c(:,k+1:m) = zero
1324 else if (beta /= one) then
1325 c(:,k+1:m) = beta * c(:,k+1:m)
1326 end if
1327 end if
1328 end if
1329
1330 ! Formatting
1331100 format(a, i0, a)
1332 end subroutine
1333
1334! ------------------------------------------------------------------------------
1335 module subroutine diag_mtx_mult_mtx2_mix(lside, alpha, a, b, err)
1336 ! Arguments
1337 logical, intent(in) :: lside
1338 complex(real64), intent(in) :: alpha
1339 real(real64), intent(in), dimension(:) :: a
1340 complex(real64), intent(inout), dimension(:,:) :: b
1341 class(errors), intent(inout), optional, target :: err
1342
1343 ! Parameters
1344 complex(real64), parameter :: zero = (0.0d0, 0.0d0)
1345 complex(real64), parameter :: one = (1.0d0, 0.0d0)
1346
1347 ! Local Variables
1348 integer(int32) :: i, m, n, k
1349 complex(real64) :: temp
1350 class(errors), pointer :: errmgr
1351 type(errors), target :: deferr
1352
1353 ! Initialization
1354 m = size(b, 1)
1355 n = size(b, 2)
1356 k = size(a)
1357 if (present(err)) then
1358 errmgr => err
1359 else
1360 errmgr => deferr
1361 end if
1362
1363 ! Input Check
1364 if ((lside .and. k > m) .or. (.not.lside .and. k > n)) then
1365 ! ERROR: One of the input arrays is not sized correctly
1366 call errmgr%report_error("diag_mtx_mult_mtx2_cmplx", &
1367 "Input number 3 is not sized correctly.", &
1368 la_array_size_error)
1369 return
1370 end if
1371
1372 ! Process
1373 if (lside) then
1374 ! Compute B = alpha * A * B
1375 do i = 1, k
1376 temp = alpha * a(i)
1377 b(i,:) = temp * b(i,:)
1378 end do
1379 if (m > k) b(k+1:m,:) = zero
1380 else
1381 ! Compute B = alpha * B * A
1382 do i = 1, k
1383 temp = alpha * a(i)
1384 b(:,i) = temp * b(:,i)
1385 end do
1386 if (n > k) b(:,k+1:n) = zero
1387 end if
1388 end subroutine
1389
1390! ------------------------------------------------------------------------------
1391 module subroutine diag_mtx_sparse_mult(lside, alpha, a, b, err)
1392 ! Arguments
1393 logical, intent(in) :: lside
1394 real(real64), intent(in) :: alpha
1395 real(real64), intent(in), dimension(:) :: a
1396 type(csr_matrix), intent(inout) :: b
1397 class(errors), intent(inout), optional, target :: err
1398
1399 ! Local Variables
1400 integer(int32) :: ii, k, k1, k2, nrow
1401 real(real64) :: scal
1402 class(errors), pointer :: errmgr
1403 type(errors), target :: deferr
1404
1405 ! Initialization
1406 nrow = size(b, 1)
1407 if (present(err)) then
1408 errmgr => err
1409 else
1410 errmgr => deferr
1411 end if
1412
1413 ! Input Check
1414 if (lside) then
1415 if (size(a) /= nrow) then
1416 call errmgr%report_error("diag_mtx_sparse_mult", &
1417 "Inner matrix dimension error.", la_array_size_error)
1418 return
1419 end if
1420 else
1421 if (size(a) /= size(b, 2)) then
1422 call errmgr%report_error("diag_mtx_sparse_mult", &
1423 "Inner matrix dimension error.", la_array_size_error)
1424 return
1425 end if
1426 end if
1427
1428 ! Process
1429 if (lside) then
1430 ! Compute B = DIAG * B
1431 do ii = 1, nrow
1432 k1 = b%row_indices(ii)
1433 k2 = b%row_indices(ii+1) - 1
1434 if (alpha == 1.0d0) then
1435 scal = a(ii)
1436 else
1437 scal = alpha * a(ii)
1438 end if
1439 do k = k1, k2
1440 b%values(k) = b%values(k) * scal
1441 end do
1442 end do
1443 else
1444 ! Compute B = B * DIAG
1445 do ii = 1, nrow
1446 k1 = b%row_indices(ii)
1447 k2 = b%row_indices(ii+1) - 1
1448 if (alpha == 1.0d0) then
1449 do k = k1, k2
1450 b%values(k) = b%values(k) * a(b%column_indices(k))
1451 end do
1452 else
1453 do k = k1, k2
1454 b%values(k) = alpha * b%values(k) * a(b%column_indices(k))
1455 end do
1456 end if
1457 end do
1458 end if
1459 end subroutine
1460
1461! ******************************************************************************
1462! BASIC OPERATION ROUTINES
1463! ------------------------------------------------------------------------------
1464 pure module function trace_dbl(x) result(y)
1465 ! Arguments
1466 real(real64), intent(in), dimension(:,:) :: x
1467 real(real64) :: y
1468
1469 ! Parameters
1470 real(real64), parameter :: zero = 0.0d0
1471
1472 ! Local Variables
1473 integer(int32) :: i, m, n, mn
1474
1475 ! Initialization
1476 y = zero
1477 m = size(x, 1)
1478 n = size(x, 2)
1479 mn = min(m, n)
1480
1481 ! Process
1482 do i = 1, mn
1483 y = y + x(i,i)
1484 end do
1485 end function
1486
1487! ------------------------------------------------------------------------------
1488 pure module function trace_cmplx(x) result(y)
1489 ! Arguments
1490 complex(real64), intent(in), dimension(:,:) :: x
1491 complex(real64) :: y
1492
1493 ! Parameters
1494 complex(real64), parameter :: zero = (0.0d0, 0.0d0)
1495
1496 ! Local Variables
1497 integer(int32) :: i, m, n, mn
1498
1499 ! Initialization
1500 y = zero
1501 m = size(x, 1)
1502 n = size(x, 2)
1503 mn = min(m, n)
1504
1505 ! Process
1506 do i = 1, mn
1507 y = y + x(i,i)
1508 end do
1509 end function
1510
1511! ------------------------------------------------------------------------------
1512 module function mtx_rank_dbl(a, tol, work, olwork, err) result(rnk)
1513 ! Arguments
1514 real(real64), intent(inout), dimension(:,:) :: a
1515 real(real64), intent(in), optional :: tol
1516 real(real64), intent(out), target, optional, dimension(:) :: work
1517 integer(int32), intent(out), optional :: olwork
1518 class(errors), intent(inout), optional, target :: err
1519 integer(int32) :: rnk
1520
1521 ! Local Variables
1522 integer(int32) :: i, m, n, mn, istat, lwork, flag
1523 real(real64), pointer, dimension(:) :: wptr, s, w
1524 real(real64), allocatable, target, dimension(:) :: wrk
1525 real(real64) :: t, tref, smlnum
1526 real(real64), dimension(1) :: dummy, temp
1527 class(errors), pointer :: errmgr
1528 type(errors), target :: deferr
1529 character(len = :), allocatable :: errmsg
1530
1531 ! Initialization
1532 m = size(a, 1)
1533 n = size(a, 2)
1534 mn = min(m, n)
1535 smlnum = dlamch('s')
1536 rnk = 0
1537 if (present(err)) then
1538 errmgr => err
1539 else
1540 errmgr => deferr
1541 end if
1542
1543 ! Workspace Query
1544 !call svd(a, a(1:mn,1), olwork = lwork)
1545 call dgesvd('N', 'N', m, n, a, m, dummy, dummy, m, dummy, n, temp, &
1546 -1, flag)
1547 lwork = int(temp(1), int32) + mn
1548 if (present(olwork)) then
1549 olwork = lwork
1550 return
1551 end if
1552
1553 ! Local Memory Allocation
1554 if (present(work)) then
1555 if (size(work) < lwork) then
1556 ! ERROR: WORK not sized correctly
1557 call errmgr%report_error("mtx_rank", &
1558 "Incorrectly sized input array WORK, argument 5.", &
1559 la_array_size_error)
1560 return
1561 end if
1562 wptr => work(1:lwork)
1563 else
1564 allocate(wrk(lwork), stat = istat)
1565 if (istat /= 0) then
1566 ! ERROR: Out of memory
1567 call errmgr%report_error("mtx_rank", &
1568 "Insufficient memory available.", &
1569 la_out_of_memory_error)
1570 return
1571 end if
1572 wptr => wrk
1573 end if
1574 s => wptr(1:mn)
1575 w => wptr(mn+1:lwork)
1576
1577 ! Compute the singular values of A
1578 call dgesvd('N', 'N', m, n, a, m, s, dummy, m, dummy, n, w, &
1579 lwork - mn, flag)
1580 if (flag > 0) then
1581 allocate(character(len = 256) :: errmsg)
1582 write(errmsg, 100) flag, " superdiagonals could not " // &
1583 "converge to zero as part of the QR iteration process."
1584 call errmgr%report_warning("mtx_rank", errmsg, la_convergence_error)
1585 end if
1586
1587 ! Determine the threshold tolerance for the singular values such that
1588 ! singular values less than the threshold result in zero when inverted.
1589 tref = max(m, n) * epsilon(t) * s(1)
1590 if (present(tol)) then
1591 t = tol
1592 else
1593 t = tref
1594 end if
1595 if (t < smlnum) then
1596 ! ! The supplied tolerance is too small, simply fall back to the
1597 ! ! default, but issue a warning to the user
1598 ! t = tref
1599 ! call report_warning("mtx_rank", "The supplied tolerance was " // &
1600 ! "smaller than a value that would result in an overflow " // &
1601 ! "condition, or is negative; therefore, the tolerance has " // &
1602 ! "been reset to its default value.")
1603 end if
1604
1605 ! Count the singular values that are larger than the tolerance value
1606 do i = 1, mn
1607 if (s(i) < t) exit
1608 rnk = rnk + 1
1609 end do
1610
1611 ! Formatting
1612100 format(i0, a)
1613 end function
1614
1615! ------------------------------------------------------------------------------
1616 module function mtx_rank_cmplx(a, tol, work, olwork, rwork, err) result(rnk)
1617 ! Arguments
1618 complex(real64), intent(inout), dimension(:,:) :: a
1619 real(real64), intent(in), optional :: tol
1620 complex(real64), intent(out), target, optional, dimension(:) :: work
1621 integer(int32), intent(out), optional :: olwork
1622 real(real64), intent(out), target, optional, dimension(:) :: rwork
1623 class(errors), intent(inout), optional, target :: err
1624 integer(int32) :: rnk
1625
1626 ! External Function Interfaces
1627 interface
1628 function dlamch(cmach) result(x)
1629 use, intrinsic :: iso_fortran_env, only : real64
1630 character, intent(in) :: cmach
1631 real(real64) :: x
1632 end function
1633 end interface
1634
1635 ! Local Variables
1636 integer(int32) :: i, m, n, mn, istat, lwork, flag, lrwork
1637 real(real64), pointer, dimension(:) :: s, rwptr, rw
1638 real(real64), allocatable, target, dimension(:) :: rwrk
1639 complex(real64), allocatable, target, dimension(:) :: wrk
1640 complex(real64), pointer, dimension(:) :: wptr
1641 real(real64) :: t, tref, smlnum
1642 real(real64), dimension(1) :: dummy
1643 complex(real64), dimension(1) :: cdummy, temp
1644 class(errors), pointer :: errmgr
1645 type(errors), target :: deferr
1646 character(len = :), allocatable :: errmsg
1647
1648 ! Initialization
1649 m = size(a, 1)
1650 n = size(a, 2)
1651 mn = min(m, n)
1652 lrwork = 6 * mn
1653 smlnum = dlamch('s')
1654 rnk = 0
1655 if (present(err)) then
1656 errmgr => err
1657 else
1658 errmgr => deferr
1659 end if
1660
1661 ! Workspace Query
1662 call zgesvd('N', 'N', m, n, a, m, dummy, cdummy, m, cdummy, n, temp, &
1663 -1, dummy, flag)
1664 lwork = int(temp(1), int32)
1665 if (present(olwork)) then
1666 olwork = lwork
1667 return
1668 end if
1669
1670 ! Local Memory Allocation
1671 if (present(work)) then
1672 if (size(work) < lwork) then
1673 ! ERROR: WORK not sized correctly
1674 call errmgr%report_error("mtx_rank_cmplx", &
1675 "Incorrectly sized input array WORK, argument 5.", &
1676 la_array_size_error)
1677 return
1678 end if
1679 wptr => work(1:lwork)
1680 else
1681 allocate(wrk(lwork), stat = istat)
1682 if (istat /= 0) then
1683 ! ERROR: Out of memory
1684 call errmgr%report_error("mtx_rank_cmplx", &
1685 "Insufficient memory available.", &
1686 la_out_of_memory_error)
1687 return
1688 end if
1689 wptr => wrk
1690 end if
1691
1692 if (present(rwork)) then
1693 if (size(rwork) < lrwork) then
1694 ! ERROR: RWORK not sized correctly
1695 call errmgr%report_error("mtx_rank_cmplx", &
1696 "Incorrectly sized input array RWORK.", &
1697 la_array_size_error)
1698 return
1699 end if
1700 rwptr => rwork(1:lrwork)
1701 else
1702 allocate(rwrk(lrwork), stat = istat)
1703 if (istat /= 0) then
1704 end if
1705 rwptr => rwrk
1706 end if
1707 s => rwptr(1:mn)
1708 rw => rwptr(mn+1:lrwork)
1709
1710 ! Compute the singular values of A
1711 call zgesvd('N', 'N', m, n, a, m, s, cdummy, m, cdummy, n, wptr, &
1712 lwork - mn, rw, flag)
1713 if (flag > 0) then
1714 allocate(character(len = 256) :: errmsg)
1715 write(errmsg, 100) flag, " superdiagonals could not " // &
1716 "converge to zero as part of the QR iteration process."
1717 call errmgr%report_warning("mtx_rank_cmplx", errmsg, la_convergence_error)
1718 end if
1719
1720 ! Determine the threshold tolerance for the singular values such that
1721 ! singular values less than the threshold result in zero when inverted.
1722 tref = max(m, n) * epsilon(t) * s(1)
1723 if (present(tol)) then
1724 t = tol
1725 else
1726 t = tref
1727 end if
1728 if (t < smlnum) then
1729 ! ! The supplied tolerance is too small, simply fall back to the
1730 ! ! default, but issue a warning to the user
1731 ! t = tref
1732 ! call report_warning("mtx_rank", "The supplied tolerance was " // &
1733 ! "smaller than a value that would result in an overflow " // &
1734 ! "condition, or is negative; therefore, the tolerance has " // &
1735 ! "been reset to its default value.")
1736 end if
1737
1738 ! Count the singular values that are larger than the tolerance value
1739 do i = 1, mn
1740 if (s(i) < t) exit
1741 rnk = rnk + 1
1742 end do
1743
1744 ! Formatting
1745100 format(i0, a)
1746 end function
1747
1748! ------------------------------------------------------------------------------
1749 module function det_dbl(a, iwork, err) result(x)
1750 ! Arguments
1751 real(real64), intent(inout), dimension(:,:) :: a
1752 integer(int32), intent(out), target, optional, dimension(:) :: iwork
1753 class(errors), intent(inout), optional, target :: err
1754 real(real64) :: x
1755
1756 ! Parameters
1757 real(real64), parameter :: zero = 0.0d0
1758 real(real64), parameter :: one = 1.0d0
1759 real(real64), parameter :: ten = 1.0d1
1760 real(real64), parameter :: p1 = 1.0d-1
1761
1762 ! Local Variables
1763 integer(int32) :: i, ep, n, istat, flag
1764 integer(int32), pointer, dimension(:) :: ipvt
1765 integer(int32), allocatable, target, dimension(:) :: iwrk
1766 real(real64) :: temp
1767 class(errors), pointer :: errmgr
1768 type(errors), target :: deferr
1769
1770 ! Initialization
1771 n = size(a, 1)
1772 x = zero
1773 if (present(err)) then
1774 errmgr => err
1775 else
1776 errmgr => deferr
1777 end if
1778
1779 ! Input Check
1780 if (size(a, 2) /= n) then
1781 call errmgr%report_error("det", &
1782 "The supplied matrix must be square.", la_array_size_error)
1783 return
1784 end if
1785
1786 ! Local Memory Allocation
1787 if (present(iwork)) then
1788 if (size(iwork) < n) then
1789 ! ERROR: WORK not sized correctly
1790 call errmgr%report_error("det", &
1791 "Incorrectly sized input array IWORK, argument 2.", &
1792 la_array_size_error)
1793 return
1794 end if
1795 ipvt => iwork(1:n)
1796 else
1797 allocate(iwrk(n), stat = istat)
1798 if (istat /= 0) then
1799 ! ERROR: Out of memory
1800 call errmgr%report_error("det", &
1801 "Insufficient memory available.", &
1802 la_out_of_memory_error)
1803 return
1804 end if
1805 ipvt => iwrk
1806 end if
1807
1808 ! Compute the LU factorization of A
1809 call dgetrf(n, n, a, n, ipvt, flag)
1810 if (flag > 0) then
1811 ! A singular matrix has a determinant of zero
1812 x = zero
1813 return
1814 end if
1815
1816 ! Compute the product of the diagonal of A
1817 temp = one
1818 ep = 0
1819 do i = 1, n
1820 if (ipvt(i) /= i) temp = -temp
1821
1822 temp = a(i,i) * temp
1823 if (temp == zero) then
1824 x = zero
1825 exit
1826 end if
1827
1828 do while (abs(temp) < one)
1829 temp = ten * temp
1830 ep = ep - 1
1831 end do
1832
1833 do while (abs(temp) > ten)
1834 temp = p1 * temp
1835 ep = ep + 1
1836 end do
1837 end do
1838 x = temp * ten**ep
1839 end function
1840
1841! ------------------------------------------------------------------------------
1842 module function det_cmplx(a, iwork, err) result(x)
1843 ! Arguments
1844 complex(real64), intent(inout), dimension(:,:) :: a
1845 integer(int32), intent(out), target, optional, dimension(:) :: iwork
1846 class(errors), intent(inout), optional, target :: err
1847 complex(real64) :: x
1848
1849 ! Parameters
1850 complex(real64), parameter :: zero = (0.0d0, 0.0d0)
1851 complex(real64), parameter :: one = (1.0d0, 0.0d0)
1852 complex(real64), parameter :: ten = (1.0d1, 0.0d0)
1853 complex(real64), parameter :: p1 = (1.0d-1, 0.0d0)
1854 real(real64), parameter :: real_one = 1.0d0
1855 real(real64), parameter :: real_ten = 1.0d1
1856
1857 ! Local Variables
1858 integer(int32) :: i, ep, n, istat, flag
1859 integer(int32), pointer, dimension(:) :: ipvt
1860 integer(int32), allocatable, target, dimension(:) :: iwrk
1861 complex(real64) :: temp
1862 class(errors), pointer :: errmgr
1863 type(errors), target :: deferr
1864
1865 ! Initialization
1866 n = size(a, 1)
1867 x = zero
1868 if (present(err)) then
1869 errmgr => err
1870 else
1871 errmgr => deferr
1872 end if
1873
1874 ! Input Check
1875 if (size(a, 2) /= n) then
1876 call errmgr%report_error("det_cmplx", &
1877 "The supplied matrix must be square.", la_array_size_error)
1878 return
1879 end if
1880
1881 ! Local Memory Allocation
1882 if (present(iwork)) then
1883 if (size(iwork) < n) then
1884 ! ERROR: WORK not sized correctly
1885 call errmgr%report_error("det_cmplx", &
1886 "Incorrectly sized input array IWORK, argument 2.", &
1887 la_array_size_error)
1888 return
1889 end if
1890 ipvt => iwork(1:n)
1891 else
1892 allocate(iwrk(n), stat = istat)
1893 if (istat /= 0) then
1894 ! ERROR: Out of memory
1895 call errmgr%report_error("det_cmplx", &
1896 "Insufficient memory available.", &
1897 la_out_of_memory_error)
1898 return
1899 end if
1900 ipvt => iwrk
1901 end if
1902
1903 ! Compute the LU factorization of A
1904 call zgetrf(n, n, a, n, ipvt, flag)
1905 if (flag > 0) then
1906 ! A singular matrix has a determinant of zero
1907 x = zero
1908 return
1909 end if
1910
1911 ! Compute the product of the diagonal of A
1912 temp = one
1913 ep = 0
1914 do i = 1, n
1915 if (ipvt(i) /= i) temp = -temp
1916
1917 temp = a(i,i) * temp
1918 if (temp == zero) then
1919 x = zero
1920 exit
1921 end if
1922
1923 do while (abs(temp) < real_one)
1924 temp = ten * temp
1925 ep = ep - 1
1926 end do
1927
1928 do while (abs(temp) > real_ten)
1929 temp = p1 * temp
1930 ep = ep + 1
1931 end do
1932 end do
1933 x = temp * ten**ep
1934 end function
1935
1936! ******************************************************************************
1937! ARRAY SWAPPING ROUTINE
1938! ------------------------------------------------------------------------------
1939 module subroutine swap_dbl(x, y, err)
1940 ! Arguments
1941 real(real64), intent(inout), dimension(:) :: x, y
1942 class(errors), intent(inout), optional, target :: err
1943
1944 ! Local Variables
1945 integer(int32) :: i, n
1946 real(real64) :: temp
1947 class(errors), pointer :: errmgr
1948 type(errors), target :: deferr
1949
1950 ! Initialization
1951 n = size(x)
1952 if (present(err)) then
1953 errmgr => err
1954 else
1955 errmgr => deferr
1956 end if
1957
1958 ! Input Check
1959 if (size(y) /= n) then
1960 call errmgr%report_error("swap", &
1961 "The input arrays are not the same size.", &
1962 la_array_size_error)
1963 return
1964 end if
1965
1966 ! Process
1967 do i = 1, n
1968 temp = x(i)
1969 x(i) = y(i)
1970 y(i) = temp
1971 end do
1972 end subroutine
1973
1974! ------------------------------------------------------------------------------
1975 module subroutine swap_cmplx(x, y, err)
1976 ! Arguments
1977 complex(real64), intent(inout), dimension(:) :: x, y
1978 class(errors), intent(inout), optional, target :: err
1979
1980 ! Local Variables
1981 integer(int32) :: i, n
1982 complex(real64) :: temp
1983 class(errors), pointer :: errmgr
1984 type(errors), target :: deferr
1985
1986 ! Initialization
1987 n = size(x)
1988 if (present(err)) then
1989 errmgr => err
1990 else
1991 errmgr => deferr
1992 end if
1993
1994 ! Input Check
1995 if (size(y) /= n) then
1996 call errmgr%report_error("swap_cmplx", &
1997 "The input arrays are not the same size.", &
1998 la_array_size_error)
1999 return
2000 end if
2001
2002 ! Process
2003 do i = 1, n
2004 temp = x(i)
2005 x(i) = y(i)
2006 y(i) = temp
2007 end do
2008 end subroutine
2009
2010! ******************************************************************************
2011! ARRAY MULTIPLICIATION ROUTINES
2012! ------------------------------------------------------------------------------
2013 module subroutine recip_mult_array_dbl(a, x)
2014 ! Arguments
2015 real(real64), intent(in) :: a
2016 real(real64), intent(inout), dimension(:) :: x
2017
2018 ! External Function Interfaces
2019 interface
2020 function dlamch(cmach) result(x)
2021 use, intrinsic :: iso_fortran_env, only : real64
2022 character, intent(in) :: cmach
2023 real(real64) :: x
2024 end function
2025 end interface
2026
2027 ! Parameters
2028 real(real64), parameter :: zero = 0.0d0
2029 real(real64), parameter :: one = 1.0d0
2030 real(real64), parameter :: twotho = 2.0d3
2031
2032 ! Local Variables
2033 logical :: done
2034 real(real64) :: bignum, cden, cden1, cnum, cnum1, mul, smlnum
2035
2036 ! Initialization
2037 smlnum = dlamch('s')
2038 bignum = one / smlnum
2039 if (log10(bignum) > twotho) then
2040 smlnum = sqrt(smlnum)
2041 bignum = sqrt(bignum)
2042 end if
2043
2044 ! Initialize the denominator to A, and the numerator to ONE
2045 cden = a
2046 cnum = one
2047
2048 ! Process
2049 do
2050 cden1 = cden * smlnum
2051 cnum1 = cnum / bignum
2052 if (abs(cden1) > abs(cnum) .and. cnum /= zero) then
2053 mul = smlnum
2054 done = .false.
2055 cden = cden1
2056 else if (abs(cnum1) > abs(cden)) then
2057 mul = bignum
2058 done = .false.
2059 cnum = cnum1
2060 else
2061 mul = cnum / cden
2062 done = .true.
2063 end if
2064
2065 ! Scale the vector X by MUL
2066 x = mul * x
2067
2068 ! Exit if done
2069 if (done) exit
2070 end do
2071 end subroutine
2072
2073! ******************************************************************************
2074! TRIANGULAR MATRIX MULTIPLICATION ROUTINES
2075! ------------------------------------------------------------------------------
2076 module subroutine tri_mtx_mult_dbl(upper, alpha, a, beta, b, err)
2077 ! Arguments
2078 logical, intent(in) :: upper
2079 real(real64), intent(in) :: alpha, beta
2080 real(real64), intent(in), dimension(:,:) :: a
2081 real(real64), intent(inout), dimension(:,:) :: b
2082 class(errors), intent(inout), optional, target :: err
2083
2084 ! Parameters
2085 real(real64), parameter :: zero = 0.0d0
2086
2087 ! Local Variables
2088 integer(int32) :: i, j, k, n, d1, d2, flag
2089 real(real64) :: temp
2090 class(errors), pointer :: errmgr
2091 type(errors), target :: deferr
2092 character(len = :), allocatable :: errmsg
2093
2094 ! Initialization
2095 n = size(a, 1)
2096 d1 = n
2097 d2 = n
2098 if (present(err)) then
2099 errmgr => err
2100 else
2101 errmgr => deferr
2102 end if
2103
2104 ! Input Check
2105 flag = 0
2106 if (size(a, 2) /= n) then
2107 flag = 3
2108 d2 = size(a, 2)
2109 else if (size(b, 1) /= n .or. size(b, 2) /= n) then
2110 flag = 5
2111 d1 = size(b, 1)
2112 d2 = size(b, 2)
2113 end if
2114 if (flag /= 0) then
2115 ! ERROR: Incorrectly sized matrix
2116 allocate(character(len = 256) :: errmsg)
2117 write(errmsg, 100) "The matrix at input ", flag, &
2118 " was not sized appropriately. A matrix of ", n, "-by-", n, &
2119 "was expected, but a matrix of ", d1, "-by-", d2, " was found."
2120 call errmgr%report_error("tri_mtx_mult_dbl", trim(errmsg), &
2121 la_array_size_error)
2122 return
2123 end if
2124
2125 ! Process
2126 if (upper) then
2127 ! Form: B = alpha * A**T * A + beta * B
2128 if (beta == zero) then
2129 do j = 1, n
2130 do i = 1, j
2131 temp = zero
2132 do k = 1, j
2133 temp = temp + a(k,i) * a(k,j)
2134 end do
2135 temp = alpha * temp
2136 b(i,j) = temp
2137 if (i /= j) b(j,i) = temp
2138 end do
2139 end do
2140 else
2141 do j = 1, n
2142 do i = 1, j
2143 temp = zero
2144 do k = 1, j
2145 temp = temp + a(k,i) * a(k,j)
2146 end do
2147 temp = alpha * temp
2148 b(i,j) = temp + beta * b(i,j)
2149 if (i /= j) b(j,i) = temp + beta * b(j,i)
2150 end do
2151 end do
2152 end if
2153 else
2154 ! Form: B = alpha * A * A**T + beta * B
2155 if (beta == zero) then
2156 do j = 1, n
2157 do i = j, n
2158 temp = zero
2159 do k = 1, j
2160 temp = temp + a(i,k) * a(j,k)
2161 end do
2162 temp = alpha * temp
2163 b(i,j) = temp
2164 if (i /= j) b(j,i) = temp
2165 end do
2166 end do
2167 else
2168 do j = 1, n
2169 do i = j, n
2170 temp = zero
2171 do k = 1, j
2172 temp = temp + a(i,k) * a(j,k)
2173 end do
2174 temp = alpha * temp
2175 b(i,j) = temp + beta * b(i,j)
2176 if (i /= j) b(j,i) = temp + beta * b(j,i)
2177 end do
2178 end do
2179 end if
2180 end if
2181
2182 ! Formatting
2183100 format(a, i0, a, i0, a, i0, a, i0, a, i0, a)
2184 end subroutine
2185
2186! ------------------------------------------------------------------------------
2187 module subroutine tri_mtx_mult_cmplx(upper, alpha, a, beta, b, err)
2188 ! Arguments
2189 logical, intent(in) :: upper
2190 complex(real64), intent(in) :: alpha, beta
2191 complex(real64), intent(in), dimension(:,:) :: a
2192 complex(real64), intent(inout), dimension(:,:) :: b
2193 class(errors), intent(inout), optional, target :: err
2194
2195 ! Parameters
2196 complex(real64), parameter :: zero = (0.0d0, 0.0d0)
2197
2198 ! Local Variables
2199 integer(int32) :: i, j, k, n, d1, d2, flag
2200 complex(real64) :: temp
2201 class(errors), pointer :: errmgr
2202 type(errors), target :: deferr
2203 character(len = :), allocatable :: errmsg
2204
2205 ! Initialization
2206 n = size(a, 1)
2207 d1 = n
2208 d2 = n
2209 if (present(err)) then
2210 errmgr => err
2211 else
2212 errmgr => deferr
2213 end if
2214
2215 ! Input Check
2216 flag = 0
2217 if (size(a, 2) /= n) then
2218 flag = 3
2219 d2 = size(a, 2)
2220 else if (size(b, 1) /= n .or. size(b, 2) /= n) then
2221 flag = 5
2222 d1 = size(b, 1)
2223 d2 = size(b, 2)
2224 end if
2225 if (flag /= 0) then
2226 ! ERROR: Incorrectly sized matrix
2227 allocate(character(len = 256) :: errmsg)
2228 write(errmsg, 100) "The matrix at input ", flag, &
2229 " was not sized appropriately. A matrix of ", n, "-by-", n, &
2230 "was expected, but a matrix of ", d1, "-by-", d2, " was found."
2231 call errmgr%report_error("tri_mtx_mult_cmplx", trim(errmsg), &
2232 la_array_size_error)
2233 return
2234 end if
2235
2236 ! Process
2237 if (upper) then
2238 ! Form: B = alpha * A**T * A + beta * B
2239 if (beta == zero) then
2240 do j = 1, n
2241 do i = 1, j
2242 temp = zero
2243 do k = 1, j
2244 temp = temp + a(k,i) * a(k,j)
2245 end do
2246 temp = alpha * temp
2247 b(i,j) = temp
2248 if (i /= j) b(j,i) = temp
2249 end do
2250 end do
2251 else
2252 do j = 1, n
2253 do i = 1, j
2254 temp = zero
2255 do k = 1, j
2256 temp = temp + a(k,i) * a(k,j)
2257 end do
2258 temp = alpha * temp
2259 b(i,j) = temp + beta * b(i,j)
2260 if (i /= j) b(j,i) = temp + beta * b(j,i)
2261 end do
2262 end do
2263 end if
2264 else
2265 ! Form: B = alpha * A * A**T + beta * B
2266 if (beta == zero) then
2267 do j = 1, n
2268 do i = j, n
2269 temp = zero
2270 do k = 1, j
2271 temp = temp + a(i,k) * a(j,k)
2272 end do
2273 temp = alpha * temp
2274 b(i,j) = temp
2275 if (i /= j) b(j,i) = temp
2276 end do
2277 end do
2278 else
2279 do j = 1, n
2280 do i = j, n
2281 temp = zero
2282 do k = 1, j
2283 temp = temp + a(i,k) * a(j,k)
2284 end do
2285 temp = alpha * temp
2286 b(i,j) = temp + beta * b(i,j)
2287 if (i /= j) b(j,i) = temp + beta * b(j,i)
2288 end do
2289 end do
2290 end if
2291 end if
2292
2293 ! Formatting
2294100 format(a, i0, a, i0, a, i0, a, i0, a, i0, a)
2295 end subroutine
2296
2297! ******************************************************************************
2298! BANDED MATRIX MULTIPLICATION ROUTINES
2299! ------------------------------------------------------------------------------
2300 module subroutine band_mtx_vec_mult_dbl(trans, kl, ku, alpha, a, x, beta, &
2301 y, err)
2302 ! Arguments
2303 logical, intent(in) :: trans
2304 integer(int32), intent(in) :: kl, ku
2305 real(real64), intent(in) :: alpha, beta
2306 real(real64), intent(in), dimension(:,:) :: a
2307 real(real64), intent(in), dimension(:) :: x
2308 real(real64), intent(inout), dimension(:) :: y
2309 class(errors), intent(inout), optional, target :: err
2310
2311 ! Local Variables
2312 integer(int32) :: m, n
2313 class(errors), pointer :: errmgr
2314 type(errors), target :: deferr
2315
2316 ! Initialization
2317 if (present(err)) then
2318 errmgr => err
2319 else
2320 errmgr => deferr
2321 end if
2322 if (trans) then
2323 m = size(x)
2324 n = size(y)
2325 else
2326 m = size(y)
2327 n = size(x)
2328 end if
2329
2330 ! Input Checking
2331 if (kl < 0) go to 10
2332 if (ku < 0) go to 20
2333 if (size(a, 1) /= kl + ku + 1) go to 30
2334 if (size(a, 2) /= n) go to 30
2335
2336 ! Process
2337 if (trans) then
2338 call dgbmv("T", m, n, kl, ku, alpha, a, size(a, 1), x, 1, beta, y, 1)
2339 else
2340 call dgbmv("N", m, n, kl, ku, alpha, a, size(a, 1), x, 1, beta, y, 1)
2341 end if
2342
2343 ! End
2344 return
2345
2346 ! KL < 0
234710 continue
2348 call errmgr%report_error("band_mtx_vec_mult_dbl", &
2349 "The number of subdiagonals must be at least 0.", &
2350 la_invalid_input_error)
2351 return
2352
2353 ! KU < 0
235420 continue
2355 call errmgr%report_error("band_mtx_vec_mult_dbl", &
2356 "The number of superdiagonals must be at least 0.", &
2357 la_invalid_input_error)
2358 return
2359
2360 ! A is incorrectly sized
236130 continue
2362 call errmgr%report_error("band_mtx_vec_mult_dbl", &
2363 "The size of matrix A is not compatible with the other vectors.", &
2364 la_array_size_error)
2365 return
2366 end subroutine
2367
2368! ------------------------------------------------------------------------------
2369 module subroutine band_mtx_vec_mult_cmplx(trans, kl, ku, alpha, a, x, &
2370 beta, y, err)
2371 ! Arguments
2372 integer(int32), intent(in) :: trans
2373 integer(int32), intent(in) :: kl, ku
2374 complex(real64), intent(in) :: alpha, beta
2375 complex(real64), intent(in), dimension(:,:) :: a
2376 complex(real64), intent(in), dimension(:) :: x
2377 complex(real64), intent(inout), dimension(:) :: y
2378 class(errors), intent(inout), optional, target :: err
2379
2380 ! Local Variables
2381 character :: op
2382 logical :: trns
2383 integer(int32) :: m, n
2384 class(errors), pointer :: errmgr
2385 type(errors), target :: deferr
2386
2387 ! Initialization
2388 if (present(err)) then
2389 errmgr => err
2390 else
2391 errmgr => deferr
2392 end if
2393 if (trans == la_transpose) then
2394 op = "T"
2395 trns = .true.
2396 else if (trans == la_hermitian_transpose) then
2397 op = "C"
2398 trns = .true.
2399 else
2400 op = "N"
2401 trns = .false.
2402 end if
2403 if (trns) then
2404 m = size(x)
2405 n = size(y)
2406 else
2407 m = size(y)
2408 n = size(x)
2409 end if
2410
2411 ! Input Checking
2412 if (kl < 0) go to 10
2413 if (ku < 0) go to 20
2414 if (size(a, 1) /= kl + ku + 1) go to 30
2415 if (size(a, 2) /= n) go to 30
2416
2417 ! Process
2418 call zgbmv(op, m, n, kl, ku, alpha, a, size(a, 1), x, 1, beta, y, 1)
2419
2420 ! End
2421 return
2422
2423 ! KL < 0
242410 continue
2425 call errmgr%report_error("band_mtx_vec_mult_cmplx", &
2426 "The number of subdiagonals must be at least 0.", &
2427 la_invalid_input_error)
2428 return
2429
2430 ! KU < 0
243120 continue
2432 call errmgr%report_error("band_mtx_vec_mult_cmplx", &
2433 "The number of superdiagonals must be at least 0.", &
2434 la_invalid_input_error)
2435 return
2436
2437 ! A is incorrectly sized
243830 continue
2439 call errmgr%report_error("band_mtx_vec_mult_cmplx", &
2440 "The size of matrix A is not compatible with the other vectors.", &
2441 la_array_size_error)
2442 return
2443 end subroutine
2444
2445! ------------------------------------------------------------------------------
2446 module subroutine band_to_full_mtx_dbl(kl, ku, b, f, err)
2447 ! Arguments
2448 integer(int32), intent(in) :: kl, ku
2449 real(real64), intent(in), dimension(:,:) :: b
2450 real(real64), intent(out), dimension(:,:) :: f
2451 class(errors), intent(inout), optional, target :: err
2452
2453 ! Parameters
2454 real(real64), parameter :: zero = 0.0d0
2455
2456 ! Local Variables
2457 class(errors), pointer :: errmgr
2458 type(errors), target :: deferr
2459 integer(int32) :: i, j, k, m, n, i1, i2
2460
2461 ! Initialization
2462 if (present(err)) then
2463 errmgr => err
2464 else
2465 errmgr => deferr
2466 end if
2467 m = size(f, 1)
2468 n = size(f, 2)
2469
2470 ! Input Check
2471 if (kl < 0) go to 10
2472 if (ku < 0) go to 20
2473 if (size(b, 2) /= n) go to 30
2474 if (size(b, 1) /= kl + ku + 1) go to 40
2475
2476 ! Process
2477 do j = 1, n
2478 k = ku + 1 - j
2479 i1 = max(1, j - ku)
2480 i2 = min(m, j + kl)
2481 do i = 1, i1 - 1
2482 f(i,j) = zero
2483 end do
2484 do i = i1, i2
2485 f(i,j) = b(k+i,j)
2486 end do
2487 do i = i2 + 1, m
2488 f(i,j) = zero
2489 end do
2490 end do
2491
2492 ! End
2493 return
2494
2495 ! KL < 0
249610 continue
2497 call errmgr%report_error("band_to_full_mtx_dbl", &
2498 "The number of subdiagonals must be at least 0.", &
2499 la_invalid_input_error)
2500 return
2501
2502 ! KU < 0
250320 continue
2504 call errmgr%report_error("band_to_full_mtx_dbl", &
2505 "The number of superdiagonals must be at least 0.", &
2506 la_invalid_input_error)
2507 return
2508
2509 ! A is incorrectly sized
251030 continue
2511 call errmgr%report_error("band_to_full_mtx_dbl", &
2512 "The number of columns in the banded matrix does not match " // &
2513 "the number of columns in the full matrix.", &
2514 la_array_size_error)
2515 return
2516
251740 continue
2518 call errmgr%report_error("band_to_full_mtx_dbl", &
2519 "The number of rows in the banded matrix does not align with " // &
2520 "the number of sub and super-diagonals specified.", &
2521 la_array_size_error)
2522 return
2523 end subroutine
2524
2525! ------------------------------------------------------------------------------
2526 module subroutine band_to_full_mtx_cmplx(kl, ku, b, f, err)
2527 ! Arguments
2528 integer(int32), intent(in) :: kl, ku
2529 complex(real64), intent(in), dimension(:,:) :: b
2530 complex(real64), intent(out), dimension(:,:) :: f
2531 class(errors), intent(inout), optional, target :: err
2532
2533 ! Parameters
2534 complex(real64), parameter :: zero = (0.0d0, 0.0d0)
2535
2536 ! Local Variables
2537 class(errors), pointer :: errmgr
2538 type(errors), target :: deferr
2539 integer(int32) :: i, j, k, m, n, i1, i2
2540
2541 ! Initialization
2542 if (present(err)) then
2543 errmgr => err
2544 else
2545 errmgr => deferr
2546 end if
2547 m = size(f, 1)
2548 n = size(f, 2)
2549
2550 ! Input Check
2551 if (kl < 0) go to 10
2552 if (ku < 0) go to 20
2553 if (size(b, 2) /= n) go to 30
2554 if (size(b, 1) /= kl + ku + 1) go to 40
2555
2556 ! Process
2557 do j = 1, n
2558 k = ku + 1 - j
2559 i1 = max(1, j - ku)
2560 i2 = min(m, j + kl)
2561 do i = 1, i1 - 1
2562 f(i,j) = zero
2563 end do
2564 do i = i1, i2
2565 f(i,j) = b(k+i,j)
2566 end do
2567 do i = i2 + 1, m
2568 f(i,j) = zero
2569 end do
2570 end do
2571
2572 ! End
2573 return
2574
2575 ! KL < 0
257610 continue
2577 call errmgr%report_error("band_to_full_mtx_cmplx", &
2578 "The number of subdiagonals must be at least 0.", &
2579 la_invalid_input_error)
2580 return
2581
2582 ! KU < 0
258320 continue
2584 call errmgr%report_error("band_to_full_mtx_cmplx", &
2585 "The number of superdiagonals must be at least 0.", &
2586 la_invalid_input_error)
2587 return
2588
2589 ! A is incorrectly sized
259030 continue
2591 call errmgr%report_error("band_to_full_mtx_cmplx", &
2592 "The number of columns in the banded matrix does not match " // &
2593 "the number of columns in the full matrix.", &
2594 la_array_size_error)
2595 return
2596
259740 continue
2598 call errmgr%report_error("band_to_full_mtx_cmplx", &
2599 "The number of rows in the banded matrix does not align with " // &
2600 "the number of sub and super-diagonals specified.", &
2601 la_array_size_error)
2602 return
2603 end subroutine
2604
2605! ------------------------------------------------------------------------------
2606 module subroutine band_diag_mtx_mult_dbl(left, m, kl, ku, alpha, a, b, err)
2607 ! Arguments
2608 logical, intent(in) :: left
2609 integer(int32), intent(in) :: m, kl, ku
2610 real(real64), intent(in) :: alpha
2611 real(real64), intent(inout), dimension(:,:) :: a
2612 real(real64), intent(in), dimension(:) :: b
2613 class(errors), intent(inout), optional, target :: err
2614
2615 ! Parameters
2616 real(real64), parameter :: one = 1.0d0
2617
2618 ! Local Variables
2619 class(errors), pointer :: errmgr
2620 type(errors), target :: deferr
2621 integer(int32) :: i, i1, i2, j, k, n
2622 real(real64) :: temp
2623
2624 ! Initialization
2625 if (present(err)) then
2626 errmgr => err
2627 else
2628 errmgr => deferr
2629 end if
2630 n = size(a, 2)
2631
2632 ! Input Checking
2633 if (kl < 0) go to 10
2634 if (ku < 0) go to 20
2635 if (left) then
2636 if (size(b) /= n) go to 30
2637 else
2638 if (size(b) < m) go to 30
2639 end if
2640
2641 ! Process
2642 if (left) then
2643 ! Compute A = A * B
2644 do j = 1, n
2645 k = ku + 1 - j
2646 i1 = max(1, j - ku) + k
2647 i2 = min(m, j + kl) + k
2648 if (alpha == one) then
2649 temp = b(j)
2650 else
2651 temp = alpha * b(j)
2652 end if
2653 do i = i1, i2
2654 a(i,j) = a(i,j) * temp
2655 end do
2656 end do
2657 else
2658 ! Compute A = B * A
2659 do j = 1, n
2660 k = ku + 1 - j
2661 i1 = max(1, j - ku)
2662 i2 = min(m, j + kl)
2663 if (alpha == 1.0d0) then
2664 do i = i1, i2
2665 a(i+k,j) = a(i+k,j) * b(i)
2666 end do
2667 else
2668 do i = i1, i2
2669 a(i+k,j) = alpha * a(i+k,j) * b(i)
2670 end do
2671 end if
2672 end do
2673 end if
2674
2675
2676 ! End
2677 return
2678
2679 ! KL < 0
268010 continue
2681 call errmgr%report_error("band_diag_mtx_mult_dbl", &
2682 "The number of subdiagonals must be at least 0.", &
2683 la_invalid_input_error)
2684 return
2685
2686 ! KU < 0
268720 continue
2688 call errmgr%report_error("band_diag_mtx_mult_dbl", &
2689 "The number of superdiagonals must be at least 0.", &
2690 la_invalid_input_error)
2691 return
2692
2693 ! B is not sized correctly
269430 continue
2695 call errmgr%report_error("band_diag_mtx_mult_dbl", &
2696 "Inner matrix dimensions do not agree.", &
2697 la_array_size_error)
2698 return
2699 end subroutine
2700
2701! ------------------------------------------------------------------------------
2702 module subroutine band_diag_mtx_mult_cmplx(left, m, kl, ku, alpha, a, b, err)
2703 ! Arguments
2704 logical, intent(in) :: left
2705 integer(int32), intent(in) :: m, kl, ku
2706 complex(real64), intent(in) :: alpha
2707 complex(real64), intent(inout), dimension(:,:) :: a
2708 complex(real64), intent(in), dimension(:) :: b
2709 class(errors), intent(inout), optional, target :: err
2710
2711 ! Parameters
2712 complex(real64), parameter :: one = (1.0d0, 0.0d0)
2713
2714 ! Local Variables
2715 class(errors), pointer :: errmgr
2716 type(errors), target :: deferr
2717 integer(int32) :: i, i1, i2, j, k, n
2718 complex(real64) :: temp
2719
2720 ! Initialization
2721 if (present(err)) then
2722 errmgr => err
2723 else
2724 errmgr => deferr
2725 end if
2726 n = size(a, 2)
2727
2728 ! Input Checking
2729 if (kl < 0) go to 10
2730 if (ku < 0) go to 20
2731 if (left) then
2732 if (size(b) /= n) go to 30
2733 else
2734 if (size(b) < m) go to 30
2735 end if
2736
2737 ! Process
2738 if (left) then
2739 ! Compute A = A * B
2740 do j = 1, n
2741 k = ku + 1 - j
2742 i1 = max(1, j - ku) + k
2743 i2 = min(m, j + kl) + k
2744 if (alpha == one) then
2745 temp = b(j)
2746 else
2747 temp = alpha * b(j)
2748 end if
2749 do i = i1, i2
2750 a(i,j) = a(i,j) * temp
2751 end do
2752 end do
2753 else
2754 ! Compute A = B * A
2755 do j = 1, n
2756 k = ku + 1 - j
2757 i1 = max(1, j - ku)
2758 i2 = min(m, j + kl)
2759 if (alpha == 1.0d0) then
2760 do i = i1, i2
2761 a(i+k,j) = a(i+k,j) * b(i)
2762 end do
2763 else
2764 do i = i1, i2
2765 a(i+k,j) = alpha * a(i+k,j) * b(i)
2766 end do
2767 end if
2768 end do
2769 end if
2770
2771
2772 ! End
2773 return
2774
2775 ! KL < 0
277610 continue
2777 call errmgr%report_error("band_diag_mtx_mult_cmplx", &
2778 "The number of subdiagonals must be at least 0.", &
2779 la_invalid_input_error)
2780 return
2781
2782 ! KU < 0
278320 continue
2784 call errmgr%report_error("band_diag_mtx_mult_cmplx", &
2785 "The number of superdiagonals must be at least 0.", &
2786 la_invalid_input_error)
2787 return
2788
2789 ! B is not sized correctly
279030 continue
2791 call errmgr%report_error("band_diag_mtx_mult_cmplx", &
2792 "Inner matrix dimensions do not agree.", &
2793 la_array_size_error)
2794 return
2795 end subroutine
2796
2797! ------------------------------------------------------------------------------
2798 module subroutine banded_to_dense_dbl(m, kl, ku, a, x, err)
2799 ! Arguments
2800 integer(int32), intent(in) :: m, kl, ku
2801 real(real64), intent(in), dimension(:,:) :: a
2802 real(real64), intent(out), dimension(:,:) :: x
2803 class(errors), intent(inout), optional, target :: err
2804
2805 ! Parameters
2806 real(real64), parameter :: zero = 0.0d0
2807
2808 ! Local Variables
2809 integer(int32) :: i, j, k, n, i1, i2
2810 class(errors), pointer :: errmgr
2811 type(errors), target :: deferr
2812
2813 ! Initialization
2814 if (present(err)) then
2815 errmgr => err
2816 else
2817 errmgr => deferr
2818 end if
2819 n = size(a, 2)
2820
2821 ! Input Checking
2822 if (kl < 0 .or. ku < 0) then
2823 call errmgr%report_error("banded_to_dense_dbl", &
2824 "The bandwidth dimensions must not be negative-valued.", &
2825 la_invalid_input_error)
2826 return
2827 end if
2828 if (size(a, 1) /= kl + ku + 1) then
2829 call errmgr%report_error("banded_to_dense_dbl", "The size of " // &
2830 "the input matrix does not match the specified bandwidth.", &
2831 la_matrix_format_error)
2832 return
2833 end if
2834 if (size(x, 1) /= m .or. size(x, 2) /= n) then
2835 call errmgr%report_error("banded_to_dense_dbl", &
2836 "The output matrix dimensions are not correct.", &
2837 la_array_size_error)
2838 return
2839 end if
2840
2841 ! Process
2842 do j = 1, n
2843 k = ku + 1 - j
2844 i1 = max(1, j - ku)
2845 i2 = min(m, j + kl)
2846 x(:i1-1,j) = zero
2847 do i = i1, i2
2848 x(i, j) = a(k + i, j)
2849 end do
2850 x(i2+1:,j) = zero
2851 end do
2852 end subroutine
2853
2854! ------------------------------------------------------------------------------
2855 module subroutine banded_to_dense_cmplx(m, kl, ku, a, x, err)
2856 ! Arguments
2857 integer(int32), intent(in) :: m, kl, ku
2858 complex(real64), intent(in), dimension(:,:) :: a
2859 complex(real64), intent(out), dimension(:,:) :: x
2860 class(errors), intent(inout), optional, target :: err
2861
2862 ! Parameters
2863 complex(real64), parameter :: zero = (0.0d0, 0.0d0)
2864
2865 ! Local Variables
2866 integer(int32) :: i, j, k, n, i1, i2
2867 class(errors), pointer :: errmgr
2868 type(errors), target :: deferr
2869
2870 ! Initialization
2871 if (present(err)) then
2872 errmgr => err
2873 else
2874 errmgr => deferr
2875 end if
2876 n = size(a, 2)
2877
2878 ! Input Checking
2879 if (kl < 0 .or. ku < 0) then
2880 call errmgr%report_error("banded_to_dense_cmplx", &
2881 "The bandwidth dimensions must not be negative-valued.", &
2882 la_invalid_input_error)
2883 return
2884 end if
2885 if (size(a, 1) /= kl + ku + 1) then
2886 call errmgr%report_error("banded_to_dense_cmplx", "The size of " // &
2887 "the input matrix does not match the specified bandwidth.", &
2888 la_matrix_format_error)
2889 return
2890 end if
2891 if (size(x, 1) /= m .or. size(x, 2) /= n) then
2892 call errmgr%report_error("banded_to_dense_cmplx", &
2893 "The output matrix dimensions are not correct.", &
2894 la_array_size_error)
2895 return
2896 end if
2897
2898 ! Process
2899 do j = 1, n
2900 k = ku + 1 - j
2901 i1 = max(1, j - ku)
2902 i2 = min(m, j + kl)
2903 x(:i1-1,j) = zero
2904 do i = i1, i2
2905 x(i, j) = a(k + i, j)
2906 end do
2907 x(i2+1:,j) = zero
2908 end do
2909 end subroutine
2910
2911! ------------------------------------------------------------------------------
2912 module subroutine dense_to_banded_dbl(a, kl, ku, x, err)
2913 ! Arguments
2914 real(real64), intent(in), dimension(:,:) :: a
2915 integer(int32), intent(in) :: kl, ku
2916 real(real64), intent(out), dimension(:,:) :: x
2917 class(errors), intent(inout), optional, target :: err
2918
2919 ! Local Variables
2920 integer(int32) :: i, j, k, m, n, mm, flag
2921 class(errors), pointer :: errmgr
2922 type(errors), target :: deferr
2923
2924 ! Initialization
2925 if (present(err)) then
2926 errmgr => err
2927 else
2928 errmgr => deferr
2929 end if
2930 m = size(a, 1)
2931 n = size(a, 2)
2932 mm = kl + ku + 1
2933
2934 ! Input Check
2935 if (kl < 0 .or. ku < 0) then
2936 call errmgr%report_error("dense_to_banded_dbl", &
2937 "The bandwidth dimensions must not be negative-valued.", &
2938 la_invalid_input_error)
2939 return
2940 end if
2941 if (size(x, 1) /= mm .or. size(x, 2) /= n) then
2942 call errmgr%report_error("dense_to_banded_dbl", &
2943 "The output matrix dimensions are not correct.", &
2944 la_array_size_error)
2945 return
2946 end if
2947
2948 ! Process
2949 do j = 1, n
2950 k = ku + 1 - j
2951 do i = max(1, j - ku), min(m, j + kl)
2952 x(k + i, j) = a(i,j)
2953 end do
2954 end do
2955 end subroutine
2956
2957! ------------------------------------------------------------------------------
2958 module subroutine dense_to_banded_cmplx(a, kl, ku, x, err)
2959 ! Arguments
2960 complex(real64), intent(in), dimension(:,:) :: a
2961 integer(int32), intent(in) :: kl, ku
2962 complex(real64), intent(out), dimension(:,:) :: x
2963 class(errors), intent(inout), optional, target :: err
2964
2965 ! Local Variables
2966 integer(int32) :: i, j, k, m, n, mm, flag
2967 class(errors), pointer :: errmgr
2968 type(errors), target :: deferr
2969
2970 ! Initialization
2971 if (present(err)) then
2972 errmgr => err
2973 else
2974 errmgr => deferr
2975 end if
2976 m = size(a, 1)
2977 n = size(a, 2)
2978 mm = kl + ku + 1
2979
2980 ! Input Check
2981 if (kl < 0 .or. ku < 0) then
2982 call errmgr%report_error("dense_to_banded_cmplx", &
2983 "The bandwidth dimensions must not be negative-valued.", &
2984 la_invalid_input_error)
2985 return
2986 end if
2987 if (size(x, 1) /= mm .or. size(x, 2) /= n) then
2988 call errmgr%report_error("dense_to_banded_cmplx", &
2989 "The output matrix dimensions are not correct.", &
2990 la_array_size_error)
2991 return
2992 end if
2993
2994 ! Process
2995 do j = 1, n
2996 k = ku + 1 - j
2997 do i = max(1, j - ku), min(m, j + kl)
2998 x(k + i, j) = a(i,j)
2999 end do
3000 end do
3001 end subroutine
3002
3003! ------------------------------------------------------------------------------
3004 module subroutine extract_diagonal_dbl(a, diag, err)
3005 ! Arguments
3006 real(real64), intent(in), dimension(:,:) :: a
3007 real(real64), intent(out), dimension(:) :: diag
3008 class(errors), intent(inout), optional, target :: err
3009
3010 ! Local Variables
3011 integer(int32) :: i, m, n, mn
3012 class(errors), pointer :: errmgr
3013 type(errors), target :: deferr
3014
3015 ! Initialization
3016 if (present(err)) then
3017 errmgr => err
3018 else
3019 errmgr => deferr
3020 end if
3021 m = size(a, 1)
3022 n = size(a, 2)
3023 mn = min(m, n)
3024
3025 ! Input Checking
3026 if (size(diag) /= mn) then
3027 call errmgr%report_error("extract_diagonal_dbl", &
3028 "The is expected to have MIN(M, N) elements.", &
3029 la_array_size_error)
3030 return
3031 end if
3032
3033 ! Process
3034 do i = 1, mn
3035 diag(i) = a(i,i)
3036 end do
3037 end subroutine
3038
3039! ------------------------------------------------------------------------------
3040 module subroutine extract_diagonal_cmplx(a, diag, err)
3041 ! Arguments
3042 complex(real64), intent(in), dimension(:,:) :: a
3043 complex(real64), intent(out), dimension(:) :: diag
3044 class(errors), intent(inout), optional, target :: err
3045
3046 ! Local Variables
3047 integer(int32) :: i, m, n, mn
3048 class(errors), pointer :: errmgr
3049 type(errors), target :: deferr
3050
3051 ! Initialization
3052 if (present(err)) then
3053 errmgr => err
3054 else
3055 errmgr => deferr
3056 end if
3057 m = size(a, 1)
3058 n = size(a, 2)
3059 mn = min(m, n)
3060
3061 ! Input Checking
3062 if (size(diag) /= mn) then
3063 call errmgr%report_error("extract_diagonal_cmplx", &
3064 "The is expected to have MIN(M, N) elements.", &
3065 la_array_size_error)
3066 return
3067 end if
3068
3069 ! Process
3070 do i = 1, mn
3071 diag(i) = a(i,i)
3072 end do
3073 end subroutine
3074
3075! ------------------------------------------------------------------------------
3076end submodule
blas
A module providing explicit interfaces to BLAS routines.
Definition blas.f90:2
linalg
Provides a set of common linear algebra routines.
Definition linalg.f90:145