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_factor.f90
1! linalg_factor.f90
2
7submodule(linalg) linalg_factor
8 use lapack
9 use qrupdate
10 implicit none
11contains
12! ******************************************************************************
13! LU FACTORIZATION
14! ------------------------------------------------------------------------------
15 module subroutine lu_factor_dbl(a, ipvt, err)
16 ! Arguments
17 real(real64), intent(inout), dimension(:,:) :: a
18 integer(int32), intent(out), dimension(:) :: ipvt
19 class(errors), intent(inout), optional, target :: err
20
21 ! Local Variables
22 integer(int32) :: m, n, mn, flag
23 class(errors), pointer :: errmgr
24 type(errors), target :: deferr
25 character(len = :), allocatable :: errmsg
26
27 ! Initialization
28 m = size(a, 1)
29 n = size(a, 2)
30 mn = min(m, n)
31 if (present(err)) then
32 errmgr => err
33 else
34 errmgr => deferr
35 end if
36
37 ! Input Check
38 flag = 0
39 if (size(ipvt) /= mn) then
40 ! ERROR: IPVT not sized correctly
41 call errmgr%report_error("lu_factor_dbl", &
42 "Incorrectly sized input array IPVT, argument 2.", &
43 la_array_size_error)
44 return
45 end if
46
47 ! Compute the LU factorization by calling the LAPACK routine DGETRF
48 call dgetrf(m, n, a, m, ipvt, flag)
49
50 ! If flag > 0, the matrix is singular. Notice, flag should not be
51 ! able to be < 0 as we've already verrified inputs prior to making the
52 ! call to LAPACK
53 if (flag > 0) then
54 ! WARNING: Singular matrix
55 allocate(character(len = 256) :: errmsg)
56 write(errmsg, 100) &
57 "Singular matrix encountered (row ", flag, ")"
58 call errmgr%report_warning("lu_factor_dbl", trim(errmsg), &
59 la_singular_matrix_error)
60 end if
61
62 ! Formatting
63100 format(a, i0, a)
64 end subroutine
65
66! ------------------------------------------------------------------------------
67 module subroutine lu_factor_cmplx(a, ipvt, err)
68 ! Arguments
69 complex(real64), intent(inout), dimension(:,:) :: a
70 integer(int32), intent(out), dimension(:) :: ipvt
71 class(errors), intent(inout), optional, target :: err
72
73 ! Local Variables
74 integer(int32) :: m, n, mn, flag
75 class(errors), pointer :: errmgr
76 type(errors), target :: deferr
77 character(len = :), allocatable :: errmsg
78
79 ! Initialization
80 m = size(a, 1)
81 n = size(a, 2)
82 mn = min(m, n)
83 if (present(err)) then
84 errmgr => err
85 else
86 errmgr => deferr
87 end if
88
89 ! Input Check
90 flag = 0
91 if (size(ipvt) /= mn) then
92 ! ERROR: IPVT not sized correctly
93 call errmgr%report_error("lu_factor_cmplx", &
94 "Incorrectly sized input array IPVT, argument 2.", &
95 la_array_size_error)
96 return
97 end if
98
99 ! Compute the LU factorization by calling the LAPACK routine ZGETRF
100 call zgetrf(m, n, a, m, ipvt, flag)
101
102 ! If flag > 0, the matrix is singular. Notice, flag should not be
103 ! able to be < 0 as we've already verrified inputs prior to making the
104 ! call to LAPACK
105 if (flag > 0) then
106 ! WARNING: Singular matrix
107 allocate(character(len = 256) :: errmsg)
108 write(errmsg, 100) &
109 "Singular matrix encountered (row ", flag, ")"
110 call errmgr%report_warning("lu_factor_cmplx", trim(errmsg), &
111 la_singular_matrix_error)
112 end if
113
114 ! Formatting
115100 format(a, i0, a)
116 end subroutine
117
118! ------------------------------------------------------------------------------
119 module subroutine form_lu_all(lu, ipvt, u, p, err)
120 ! Arguments
121 real(real64), intent(inout), dimension(:,:) :: lu
122 integer(int32), intent(in), dimension(:) :: ipvt
123 real(real64), intent(out), dimension(:,:) :: u, p
124 class(errors), intent(inout), optional, target :: err
125
126 ! Local Variables
127 integer(int32) :: j, jp, n, flag
128 class(errors), pointer :: errmgr
129 type(errors), target :: deferr
130 character(len = :), allocatable :: errmsg
131
132 ! Parameters
133 real(real64), parameter :: zero = 0.0d0
134 real(real64), parameter :: one = 1.0d0
135
136 ! Initialization
137 n = size(lu, 1)
138 if (present(err)) then
139 errmgr => err
140 else
141 errmgr => deferr
142 end if
143
144 ! Input Check
145 flag = 0
146 if (size(lu, 2) /= n) then
147 flag = 1
148 else if (size(ipvt) /= n) then
149 flag = 2
150 else if (size(u, 1) /= n .or. size(u, 2) /= n) then
151 flag = 3
152 else if (size(p, 1) /= n .or. size(p, 2) /= n) then
153 flag = 4
154 end if
155 if (flag /= 0) then
156 ! One of the input arrays is not sized correctly
157 allocate(character(len = 256) :: errmsg)
158 write(errmsg, 100) "Input number ", flag, &
159 " is not sized correctly."
160 call errmgr%report_error("form_lu_all", trim(errmsg), &
161 la_array_size_error)
162 return
163 end if
164
165 ! Ensure P starts off as an identity matrix
166 call dlaset('A', n, n, zero, one, p, n)
167
168 ! Process
169 do j = 1, n
170 ! Define the pivot matrix
171 jp = ipvt(j)
172 if (j /= jp) call swap(p(j,1:n), p(jp,1:n))
173
174 ! Build L and U
175 u(1:j,j) = lu(1:j,j)
176 u(j+1:n,j) = zero
177
178 if (j > 1) lu(1:j-1,j) = zero
179 lu(j,j) = one
180 end do
181
182 ! Formatting
183100 format(a, i0, a)
184 end subroutine
185
186! ------------------------------------------------------------------------------
187 module subroutine form_lu_all_cmplx(lu, ipvt, u, p, err)
188 ! Arguments
189 complex(real64), intent(inout), dimension(:,:) :: lu
190 integer(int32), intent(in), dimension(:) :: ipvt
191 complex(real64), intent(out), dimension(:,:) :: u
192 real(real64), intent(out), dimension(:,:) :: p
193 class(errors), intent(inout), optional, target :: err
194
195 ! Local Variables
196 integer(int32) :: j, jp, n, flag
197 class(errors), pointer :: errmgr
198 type(errors), target :: deferr
199 character(len = :), allocatable :: errmsg
200
201 ! Parameters
202 real(real64), parameter :: zero = 0.0d0
203 real(real64), parameter :: one = 1.0d0
204 complex(real64), parameter :: c_zero = (0.0d0, 0.0d0)
205 complex(real64), parameter :: c_one = (1.0d0, 0.0d0)
206
207 ! Initialization
208 n = size(lu, 1)
209 if (present(err)) then
210 errmgr => err
211 else
212 errmgr => deferr
213 end if
214
215 ! Input Check
216 flag = 0
217 if (size(lu, 2) /= n) then
218 flag = 1
219 else if (size(ipvt) /= n) then
220 flag = 2
221 else if (size(u, 1) /= n .or. size(u, 2) /= n) then
222 flag = 3
223 else if (size(p, 1) /= n .or. size(p, 2) /= n) then
224 flag = 4
225 end if
226 if (flag /= 0) then
227 ! One of the input arrays is not sized correctly
228 allocate(character(len = 256) :: errmsg)
229 write(errmsg, 100) "Input number ", flag, &
230 " is not sized correctly."
231 call errmgr%report_error("form_lu_all_cmplx", trim(errmsg), &
232 la_array_size_error)
233 return
234 end if
235
236 ! Ensure P starts off as an identity matrix
237 call dlaset('A', n, n, zero, one, p, n)
238
239 ! Process
240 do j = 1, n
241 ! Define the pivot matrix
242 jp = ipvt(j)
243 if (j /= jp) call swap(p(j,1:n), p(jp,1:n))
244
245 ! Build L and U
246 u(1:j,j) = lu(1:j,j)
247 u(j+1:n,j) = c_zero
248
249 if (j > 1) lu(1:j-1,j) = c_zero
250 lu(j,j) = c_one
251 end do
252
253 ! Formatting
254100 format(a, i0, a)
255 end subroutine
256
257! ------------------------------------------------------------------------------
258 module subroutine form_lu_only(lu, u, err)
259 ! Arguments
260 real(real64), intent(inout), dimension(:,:) :: lu
261 real(real64), intent(out), dimension(:,:) :: u
262 class(errors), intent(inout), optional, target :: err
263
264 ! Local Variables
265 integer(int32) :: j, n, flag
266 class(errors), pointer :: errmgr
267 type(errors), target :: deferr
268 character(len = :), allocatable :: errmsg
269
270 ! Parameters
271 real(real64), parameter :: zero = 0.0d0
272 real(real64), parameter :: one = 1.0d0
273
274 ! Initialization
275 n = size(lu, 1)
276 if (present(err)) then
277 errmgr => err
278 else
279 errmgr => deferr
280 end if
281
282 ! Input Check
283 flag = 0
284 if (size(lu, 2) /= n) then
285 flag = 2
286 else if (size(u, 1) /= n .or. size(u, 2) /= n) then
287 flag = 3
288 end if
289 if (flag /= 0) then
290 ! One of the input arrays is not sized correctly
291 allocate(character(len = 256) :: errmsg)
292 write(errmsg, 100) "Input number ", flag, &
293 " is not sized correctly."
294 call errmgr%report_error("form_lu_only", trim(errmsg), &
295 la_array_size_error)
296 return
297 end if
298
299 ! Process
300 do j = 1, n
301 ! Build L and U
302 u(1:j,j) = lu(1:j,j)
303 u(j+1:n,j) = zero
304
305 if (j > 1) lu(1:j-1,j) = zero
306 lu(j,j) = one
307 end do
308
309 ! Formatting
310100 format(a, i0, a)
311 end subroutine
312
313! ------------------------------------------------------------------------------
314 module subroutine form_lu_only_cmplx(lu, u, err)
315 ! Arguments
316 complex(real64), intent(inout), dimension(:,:) :: lu
317 complex(real64), intent(out), dimension(:,:) :: u
318 class(errors), intent(inout), optional, target :: err
319
320 ! Local Variables
321 integer(int32) :: j, n, flag
322 class(errors), pointer :: errmgr
323 type(errors), target :: deferr
324 character(len = :), allocatable :: errmsg
325
326 ! Parameters
327 complex(real64), parameter :: zero = (0.0d0, 0.0d0)
328 complex(real64), parameter :: one = (1.0d0, 0.0d0)
329
330 ! Initialization
331 n = size(lu, 1)
332 if (present(err)) then
333 errmgr => err
334 else
335 errmgr => deferr
336 end if
337
338 ! Input Check
339 flag = 0
340 if (size(lu, 2) /= n) then
341 flag = 2
342 else if (size(u, 1) /= n .or. size(u, 2) /= n) then
343 flag = 3
344 end if
345 if (flag /= 0) then
346 ! One of the input arrays is not sized correctly
347 allocate(character(len = 256) :: errmsg)
348 write(errmsg, 100) "Input number ", flag, &
349 " is not sized correctly."
350 call errmgr%report_error("form_lu_only_cmplx", trim(errmsg), &
351 la_array_size_error)
352 return
353 end if
354
355 ! Process
356 do j = 1, n
357 ! Build L and U
358 u(1:j,j) = lu(1:j,j)
359 u(j+1:n,j) = zero
360
361 if (j > 1) lu(1:j-1,j) = zero
362 lu(j,j) = one
363 end do
364
365 ! Formatting
366100 format(a, i0, a)
367 end subroutine
368
369! ******************************************************************************
370! QR FACTORIZATION
371! ------------------------------------------------------------------------------
372 module subroutine qr_factor_no_pivot(a, tau, work, olwork, err)
373 ! Arguments
374 real(real64), intent(inout), dimension(:,:) :: a
375 real(real64), intent(out), dimension(:) :: tau
376 real(real64), intent(out), target, dimension(:), optional :: work
377 integer(int32), intent(out), optional :: olwork
378 class(errors), intent(inout), optional, target :: err
379
380 ! Local Variables
381 integer(int32) :: m, n, mn, istat, lwork, flag
382 real(real64), dimension(1) :: temp
383 real(real64), pointer, dimension(:) :: wptr
384 real(real64), allocatable, target, dimension(:) :: wrk
385 class(errors), pointer :: errmgr
386 type(errors), target :: deferr
387
388 ! Initialization
389 m = size(a, 1)
390 n = size(a, 2)
391 mn = min(m, n)
392 if (present(err)) then
393 errmgr => err
394 else
395 errmgr => deferr
396 end if
397
398 ! Input Check
399 if (size(tau) /= mn) then
400 ! ERROR: TAU not sized correctly
401 call errmgr%report_error("qr_factor_no_pivot", &
402 "Incorrectly sized input array TAU, argument 2.", &
403 la_array_size_error)
404 return
405 end if
406
407 ! Workspace Query
408 call dgeqrf(m, n, a, m, tau, temp, -1, flag)
409 lwork = int(temp(1), int32)
410 if (present(olwork)) then
411 olwork = lwork
412 return
413 end if
414
415 ! Local Memory Allocation
416 if (present(work)) then
417 if (size(work) < lwork) then
418 ! ERROR: WORK not sized correctly
419 call errmgr%report_error("qr_factor_no_pivot", &
420 "Incorrectly sized input array WORK, argument 3.", &
421 la_array_size_error)
422 return
423 end if
424 wptr => work(1:lwork)
425 else
426 allocate(wrk(lwork), stat = istat)
427 if (istat /= 0) then
428 ! ERROR: Out of memory
429 call errmgr%report_error("qr_factor_no_pivot", &
430 "Insufficient memory available.", &
431 la_out_of_memory_error)
432 return
433 end if
434 wptr => wrk
435 end if
436
437 ! Call DGEQRF
438 call dgeqrf(m, n, a, m, tau, wptr, lwork, flag)
439 end subroutine
440
441! ------------------------------------------------------------------------------
442 module subroutine qr_factor_no_pivot_cmplx(a, tau, work, olwork, err)
443 ! Arguments
444 complex(real64), intent(inout), dimension(:,:) :: a
445 complex(real64), intent(out), dimension(:) :: tau
446 complex(real64), intent(out), target, dimension(:), optional :: work
447 integer(int32), intent(out), optional :: olwork
448 class(errors), intent(inout), optional, target :: err
449
450 ! Local Variables
451 integer(int32) :: m, n, mn, istat, lwork, flag
452 complex(real64), dimension(1) :: temp
453 complex(real64), pointer, dimension(:) :: wptr
454 complex(real64), allocatable, target, dimension(:) :: wrk
455 class(errors), pointer :: errmgr
456 type(errors), target :: deferr
457
458 ! Initialization
459 m = size(a, 1)
460 n = size(a, 2)
461 mn = min(m, n)
462 if (present(err)) then
463 errmgr => err
464 else
465 errmgr => deferr
466 end if
467
468 ! Input Check
469 if (size(tau) /= mn) then
470 ! ERROR: TAU not sized correctly
471 call errmgr%report_error("qr_factor_no_pivot_cmplx", &
472 "Incorrectly sized input array TAU, argument 2.", &
473 la_array_size_error)
474 return
475 end if
476
477 ! Workspace Query
478 call zgeqrf(m, n, a, m, tau, temp, -1, flag)
479 lwork = int(temp(1), int32)
480 if (present(olwork)) then
481 olwork = lwork
482 return
483 end if
484
485 ! Local Memory Allocation
486 if (present(work)) then
487 if (size(work) < lwork) then
488 ! ERROR: WORK not sized correctly
489 call errmgr%report_error("qr_factor_no_pivot_cmplx", &
490 "Incorrectly sized input array WORK, argument 3.", &
491 la_array_size_error)
492 return
493 end if
494 wptr => work(1:lwork)
495 else
496 allocate(wrk(lwork), stat = istat)
497 if (istat /= 0) then
498 ! ERROR: Out of memory
499 call errmgr%report_error("qr_factor_no_pivot_cmplx", &
500 "Insufficient memory available.", &
501 la_out_of_memory_error)
502 return
503 end if
504 wptr => wrk
505 end if
506
507 ! Call ZGEQRF
508 call zgeqrf(m, n, a, m, tau, wptr, lwork, flag)
509 end subroutine
510
511! ------------------------------------------------------------------------------
512 module subroutine qr_factor_pivot(a, tau, jpvt, work, olwork, err)
513 ! Arguments
514 real(real64), intent(inout), dimension(:,:) :: a
515 real(real64), intent(out), dimension(:) :: tau
516 integer(int32), intent(inout), dimension(:) :: jpvt
517 real(real64), intent(out), target, dimension(:), optional :: work
518 integer(int32), intent(out), optional :: olwork
519 class(errors), intent(inout), optional, target :: err
520
521 ! Local Variables
522 integer(int32) :: m, n, mn, istat, lwork, flag
523 real(real64), dimension(1) :: temp
524 real(real64), pointer, dimension(:) :: wptr
525 real(real64), allocatable, target, dimension(:) :: wrk
526 class(errors), pointer :: errmgr
527 type(errors), target :: deferr
528 character(len = :), allocatable :: errmsg
529
530 ! Initialization
531 m = size(a, 1)
532 n = size(a, 2)
533 mn = min(m, n)
534 if (present(err)) then
535 errmgr => err
536 else
537 errmgr => deferr
538 end if
539
540 ! Input Check
541 flag = 0
542 if (size(tau) /= mn) then
543 flag = 2
544 else if (size(jpvt) /= n) then
545 flag = 3
546 end if
547 if (flag /= 0) then
548 ! ERROR: One of the input arrays is not sized correctly
549 allocate(character(len = 256) :: errmsg)
550 write(errmsg, 100) "Input number ", flag, &
551 " is not sized correctly."
552 call errmgr%report_error("qr_factor_pivot", trim(errmsg), &
553 la_array_size_error)
554 return
555 end if
556
557 ! Workspace Query
558 call dgeqp3(m, n, a, m, jpvt, tau, temp, -1, flag)
559 lwork = int(temp(1), int32)
560 if (present(olwork)) then
561 olwork = lwork
562 return
563 end if
564
565 ! Local Memory Allocation
566 if (present(work)) then
567 if (size(work) < lwork) then
568 ! ERROR: WORK not sized correctly
569 call errmgr%report_error("qr_factor_pivot", &
570 "Incorrectly sized input array WORK, argument 4.", &
571 la_array_size_error)
572 return
573 end if
574 wptr => work(1:lwork)
575 else
576 allocate(wrk(lwork), stat = istat)
577 if (istat /= 0) then
578 ! ERROR: Out of memory
579 call errmgr%report_error("qr_factor_pivot", &
580 "Insufficient memory available.", &
581 la_out_of_memory_error)
582 return
583 end if
584 wptr => wrk
585 end if
586
587 ! Call DGEQP3
588 call dgeqp3(m, n, a, m, jpvt, tau, wptr, lwork, flag)
589
590 ! End
591 if (allocated(wrk)) deallocate(wrk)
592
593 ! Formatting
594100 format(a, i0, a)
595 end subroutine
596
597! ------------------------------------------------------------------------------
598 module subroutine qr_factor_pivot_cmplx(a, tau, jpvt, work, olwork, rwork, &
599 err)
600 ! Arguments
601 complex(real64), intent(inout), dimension(:,:) :: a
602 complex(real64), intent(out), dimension(:) :: tau
603 integer(int32), intent(inout), dimension(:) :: jpvt
604 complex(real64), intent(out), target, dimension(:), optional :: work
605 integer(int32), intent(out), optional :: olwork
606 real(real64), intent(out), target, dimension(:), optional :: rwork
607 class(errors), intent(inout), optional, target :: err
608
609 ! Local Variables
610 integer(int32) :: m, n, mn, istat, lwork, flag
611 complex(real64), dimension(1) :: temp
612 complex(real64), pointer, dimension(:) :: wptr
613 complex(real64), allocatable, target, dimension(:) :: wrk
614 real(real64), pointer, dimension(:) :: rptr
615 real(real64), allocatable, target, dimension(:) :: rwrk
616 class(errors), pointer :: errmgr
617 type(errors), target :: deferr
618 character(len = :), allocatable :: errmsg
619
620 ! Initialization
621 m = size(a, 1)
622 n = size(a, 2)
623 mn = min(m, n)
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 (size(tau) /= mn) then
633 flag = 2
634 else if (size(jpvt) /= n) then
635 flag = 3
636 end if
637 if (flag /= 0) then
638 ! ERROR: One of the input arrays is not sized correctly
639 allocate(character(len = 256) :: errmsg)
640 write(errmsg, 100) "Input number ", flag, &
641 " is not sized correctly."
642 call errmgr%report_error("qr_factor_pivot_cmplx", trim(errmsg), &
643 la_array_size_error)
644 return
645 end if
646 if (present(rwork)) then
647 if (size(rwork) < 2 * n) then
648 call errmgr%report_error("qr_factor_pivot_cmplx", &
649 "Incorrectly sized input array RWORK, argument 6.", &
650 la_array_size_error)
651 return
652 end if
653 rptr => rwork(1:2*n)
654 else
655 allocate(rwrk(2 * n), stat = flag)
656 if (flag /= 0) then
657 call errmgr%report_error("qr_factor_pivot_cmplx", &
658 "Insufficient memory available.", &
659 la_out_of_memory_error)
660 return
661 end if
662 rptr => rwrk
663 end if
664
665 ! Workspace Query
666 call zgeqp3(m, n, a, m, jpvt, tau, temp, -1, rptr, flag)
667 lwork = int(temp(1), int32)
668 if (present(olwork)) then
669 olwork = lwork
670 return
671 end if
672
673 ! Local Memory Allocation
674 if (present(work)) then
675 if (size(work) < lwork) then
676 ! ERROR: WORK not sized correctly
677 call errmgr%report_error("qr_factor_pivot_cmplx", &
678 "Incorrectly sized input array WORK, argument 4.", &
679 la_array_size_error)
680 return
681 end if
682 wptr => work(1:lwork)
683 else
684 allocate(wrk(lwork), stat = istat)
685 if (istat /= 0) then
686 ! ERROR: Out of memory
687 call errmgr%report_error("qr_factor_pivot_cmplx", &
688 "Insufficient memory available.", &
689 la_out_of_memory_error)
690 return
691 end if
692 wptr => wrk
693 end if
694
695 ! Call ZGEQP3
696 call zgeqp3(m, n, a, m, jpvt, tau, wptr, lwork, rptr, flag)
697
698 ! End
699 if (allocated(wrk)) deallocate(wrk)
700
701 ! Formatting
702100 format(a, i0, a)
703 end subroutine
704
705! ------------------------------------------------------------------------------
706 module subroutine form_qr_no_pivot(r, tau, q, work, olwork, err)
707 ! Arguments
708 real(real64), intent(inout), dimension(:,:) :: r
709 real(real64), intent(in), dimension(:) :: tau
710 real(real64), intent(out), dimension(:,:) :: q
711 real(real64), intent(out), target, dimension(:), optional :: work
712 integer(int32), intent(out), optional :: olwork
713 class(errors), intent(inout), optional, target :: err
714
715 ! Parameters
716 real(real64), parameter :: zero = 0.0d0
717
718 ! Local Variables
719 integer(int32) :: j, m, n, mn, qcol, istat, flag, lwork
720 real(real64), pointer, dimension(:) :: wptr
721 real(real64), allocatable, target, dimension(:) :: wrk
722 real(real64), dimension(1) :: temp
723 class(errors), pointer :: errmgr
724 type(errors), target :: deferr
725 character(len = :), allocatable :: errmsg
726
727 ! Initialization
728 m = size(r, 1)
729 n = size(r, 2)
730 mn = min(m, n)
731 qcol = size(q, 2)
732 if (present(err)) then
733 errmgr => err
734 else
735 errmgr => deferr
736 end if
737
738 ! Input Check
739 flag = 0
740 if (size(tau) /= mn) then
741 flag = 2
742 else if (size(q, 1) /= m .or. (qcol /= m .and. qcol /= n)) then
743 flag = 3
744 else if (qcol == n .and. m < n) then
745 flag = 3
746 end if
747 if (flag /= 0) then
748 ! ERROR: One of the input arrays is not sized correctly
749 allocate(character(len = 256) :: errmsg)
750 write(errmsg, 100) "Input number ", flag, &
751 " is not sized correctly."
752 call errmgr%report_error("form_qr_no_pivot", trim(errmsg), &
753 la_array_size_error)
754 return
755 end if
756
757 ! Workspace Query
758 call dorgqr(m, qcol, mn, q, m, tau, temp, -1, flag)
759 lwork = int(temp(1), int32)
760 if (present(olwork)) then
761 olwork = lwork
762 return
763 end if
764
765 ! Local Memory Allocation
766 if (present(work)) then
767 if (size(work) < lwork) then
768 ! ERROR: WORK not sized correctly
769 call errmgr%report_error("form_qr_no_pivot", &
770 "Incorrectly sized input array WORK, argument 4.", &
771 la_array_size_error)
772 return
773 end if
774 wptr => work(1:lwork)
775 else
776 allocate(wrk(lwork), stat = istat)
777 if (istat /= 0) then
778 ! ERROR: Out of memory
779 call errmgr%report_error("form_qr_no_pivot", &
780 "Insufficient memory available.", &
781 la_out_of_memory_error)
782 return
783 end if
784 wptr => wrk
785 end if
786
787 ! Copy the sub-diagonal portion of R to Q, and then zero out the
788 ! sub-diagonal portion of R
789 do j = 1, mn
790 q(j+1:m,j) = r(j+1:m,j)
791 r(j+1:m,j) = zero
792 end do
793
794 ! Build Q - Build M-by-M or M-by-N, but M-by-N only for M >= N
795 call dorgqr(m, qcol, mn, q, m, tau, wptr, lwork, flag)
796
797 ! End
798 if (allocated(wrk)) deallocate(wrk)
799
800 ! Formatting
801100 format(a, i0, a)
802 end subroutine
803
804! ------------------------------------------------------------------------------
805 module subroutine form_qr_no_pivot_cmplx(r, tau, q, work, olwork, err)
806 ! Arguments
807 complex(real64), intent(inout), dimension(:,:) :: r
808 complex(real64), intent(in), dimension(:) :: tau
809 complex(real64), intent(out), dimension(:,:) :: q
810 complex(real64), intent(out), target, dimension(:), optional :: work
811 integer(int32), intent(out), optional :: olwork
812 class(errors), intent(inout), optional, target :: err
813
814 ! Parameters
815 complex(real64), parameter :: zero = (0.0d0, 0.0d0)
816
817 ! Local Variables
818 integer(int32) :: j, m, n, mn, qcol, istat, flag, lwork
819 complex(real64), pointer, dimension(:) :: wptr
820 complex(real64), allocatable, target, dimension(:) :: wrk
821 complex(real64), dimension(1) :: temp
822 class(errors), pointer :: errmgr
823 type(errors), target :: deferr
824 character(len = :), allocatable :: errmsg
825
826 ! Initialization
827 m = size(r, 1)
828 n = size(r, 2)
829 mn = min(m, n)
830 qcol = size(q, 2)
831 if (present(err)) then
832 errmgr => err
833 else
834 errmgr => deferr
835 end if
836
837 ! Input Check
838 flag = 0
839 if (size(tau) /= mn) then
840 flag = 2
841 else if (size(q, 1) /= m .or. (qcol /= m .and. qcol /= n)) then
842 flag = 3
843 else if (qcol == n .and. m < n) then
844 flag = 3
845 end if
846 if (flag /= 0) then
847 ! ERROR: One of the input arrays is not sized correctly
848 allocate(character(len = 256) :: errmsg)
849 write(errmsg, 100) "Input number ", flag, &
850 " is not sized correctly."
851 call errmgr%report_error("form_qr_no_pivot_cmplx", trim(errmsg), &
852 la_array_size_error)
853 return
854 end if
855
856 ! Workspace Query
857 call zungqr(m, qcol, mn, q, m, tau, temp, -1, flag)
858 lwork = int(temp(1), int32)
859 if (present(olwork)) then
860 olwork = lwork
861 return
862 end if
863
864 ! Local Memory Allocation
865 if (present(work)) then
866 if (size(work) < lwork) then
867 ! ERROR: WORK not sized correctly
868 call errmgr%report_error("form_qr_no_pivot_cmplx", &
869 "Incorrectly sized input array WORK, argument 4.", &
870 la_array_size_error)
871 return
872 end if
873 wptr => work(1:lwork)
874 else
875 allocate(wrk(lwork), stat = istat)
876 if (istat /= 0) then
877 ! ERROR: Out of memory
878 call errmgr%report_error("form_qr_no_pivot_cmplx", &
879 "Insufficient memory available.", &
880 la_out_of_memory_error)
881 return
882 end if
883 wptr => wrk
884 end if
885
886 ! Copy the sub-diagonal portion of R to Q, and then zero out the
887 ! sub-diagonal portion of R
888 do j = 1, mn
889 q(j+1:m,j) = r(j+1:m,j)
890 r(j+1:m,j) = zero
891 end do
892
893 ! Build Q - Build M-by-M or M-by-N, but M-by-N only for M >= N
894 call zungqr(m, qcol, mn, q, m, tau, wptr, lwork, flag)
895
896 ! End
897 if (allocated(wrk)) deallocate(wrk)
898
899 ! Formatting
900100 format(a, i0, a)
901 end subroutine
902
903! ------------------------------------------------------------------------------
904 module subroutine form_qr_pivot(r, tau, pvt, q, p, work, olwork, err)
905 ! Arguments
906 real(real64), intent(inout), dimension(:,:) :: r
907 real(real64), intent(in), dimension(:) :: tau
908 integer(int32), intent(in), dimension(:) :: pvt
909 real(real64), intent(out), dimension(:,:) :: q, p
910 real(real64), intent(out), target, dimension(:), optional :: work
911 integer(int32), intent(out), optional :: olwork
912 class(errors), intent(inout), optional, target :: err
913
914 ! Parameters
915 real(real64), parameter :: zero = 0.0d0
916 real(real64), parameter :: one = 1.0d0
917
918 ! Local Variables
919 integer(int32) :: j, jp, m, n, mn, flag
920 class(errors), pointer :: errmgr
921 type(errors), target :: deferr
922 character(len = :), allocatable :: errmsg
923
924 ! Initialization
925 m = size(r, 1)
926 n = size(r, 2)
927 mn = min(m, n)
928 if (present(err)) then
929 errmgr => err
930 else
931 errmgr => deferr
932 end if
933
934 ! Input Check
935 flag = 0
936 if (size(tau) /= mn) then
937 flag = 2
938 else if (size(pvt) /= n) then
939 flag = 3
940 else if (size(q, 1) /= m .or. &
941 (size(q, 2) /= m .and. size(q, 2) /= n)) then
942 flag = 4
943 else if (size(q, 2) == n .and. m < n) then
944 flag = 4
945 else if (size(p, 1) /= n .or. size(p, 2) /= n) then
946 flag = 5
947 end if
948 if (flag /= 0) then
949 ! ERROR: One of the input arrays is not sized correctly
950 allocate(character(len = 256) :: errmsg)
951 write(errmsg, 100) "Input number ", flag, &
952 " is not sized correctly."
953 call errmgr%report_error("form_qr_pivot", trim(errmsg), &
954 la_array_size_error)
955 return
956 end if
957
958 ! Generate Q and R
959 call form_qr_no_pivot(r, tau, q, work = work, olwork = olwork, &
960 err = errmgr)
961 if (present(olwork)) return ! Just a workspace query
962 if (errmgr%has_error_occurred()) return
963
964 ! Form P
965 do j = 1, n
966 jp = pvt(j)
967 p(:,j) = zero
968 p(jp,j) = one
969 end do
970
971 ! Formatting
972100 format(a, i0, a)
973 end subroutine
974
975! ------------------------------------------------------------------------------
976 module subroutine form_qr_pivot_cmplx(r, tau, pvt, q, p, work, olwork, err)
977 ! Arguments
978 complex(real64), intent(inout), dimension(:,:) :: r
979 complex(real64), intent(in), dimension(:) :: tau
980 integer(int32), intent(in), dimension(:) :: pvt
981 complex(real64), intent(out), dimension(:,:) :: q, p
982 complex(real64), intent(out), target, dimension(:), optional :: work
983 integer(int32), intent(out), optional :: olwork
984 class(errors), intent(inout), optional, target :: err
985
986 ! Parameters
987 complex(real64), parameter :: zero = (0.0d0, 0.0d0)
988 complex(real64), parameter :: one = (1.0d0, 0.0d0)
989
990 ! Local Variables
991 integer(int32) :: j, jp, m, n, mn, flag
992 class(errors), pointer :: errmgr
993 type(errors), target :: deferr
994 character(len = :), allocatable :: errmsg
995
996 ! Initialization
997 m = size(r, 1)
998 n = size(r, 2)
999 mn = min(m, n)
1000 if (present(err)) then
1001 errmgr => err
1002 else
1003 errmgr => deferr
1004 end if
1005
1006 ! Input Check
1007 flag = 0
1008 if (size(tau) /= mn) then
1009 flag = 2
1010 else if (size(pvt) /= n) then
1011 flag = 3
1012 else if (size(q, 1) /= m .or. &
1013 (size(q, 2) /= m .and. size(q, 2) /= n)) then
1014 flag = 4
1015 else if (size(q, 2) == n .and. m < n) then
1016 flag = 4
1017 else if (size(p, 1) /= n .or. size(p, 2) /= n) then
1018 flag = 5
1019 end if
1020 if (flag /= 0) then
1021 ! ERROR: One of the input arrays is not sized correctly
1022 allocate(character(len = 256) :: errmsg)
1023 write(errmsg, 100) "Input number ", flag, &
1024 " is not sized correctly."
1025 call errmgr%report_error("form_qr_pivot_cmplx", trim(errmsg), &
1026 la_array_size_error)
1027 return
1028 end if
1029
1030 ! Generate Q and R
1031 call form_qr_no_pivot_cmplx(r, tau, q, work = work, olwork = olwork, &
1032 err = errmgr)
1033 if (present(olwork)) return ! Just a workspace query
1034 if (errmgr%has_error_occurred()) return
1035
1036 ! Form P
1037 do j = 1, n
1038 jp = pvt(j)
1039 p(:,j) = zero
1040 p(jp,j) = one
1041 end do
1042
1043 ! Formatting
1044100 format(a, i0, a)
1045 end subroutine
1046
1047! ------------------------------------------------------------------------------
1048 module subroutine mult_qr_mtx(lside, trans, a, tau, c, work, olwork, err)
1049 ! Arguments
1050 logical, intent(in) :: lside, trans
1051 real(real64), intent(in), dimension(:) :: tau
1052 real(real64), intent(inout), dimension(:,:) :: a, c
1053 real(real64), intent(out), target, dimension(:), optional :: work
1054 integer(int32), intent(out), optional :: olwork
1055 class(errors), intent(inout), optional, target :: err
1056
1057 ! Parameters
1058 real(real64), parameter :: one = 1.0d0
1059
1060 ! Local Variables
1061 character :: side, t
1062 integer(int32) :: m, n, k, nrowa, istat, flag, lwork
1063 real(real64), pointer, dimension(:) :: wptr
1064 real(real64), allocatable, target, dimension(:) :: wrk
1065 real(real64), dimension(1) :: temp
1066 class(errors), pointer :: errmgr
1067 type(errors), target :: deferr
1068 character(len = :), allocatable :: errmsg
1069
1070 ! Initialization
1071 m = size(c, 1)
1072 n = size(c, 2)
1073 k = size(tau)
1074 if (lside) then
1075 side = 'L'
1076 nrowa = m
1077 else
1078 side = 'R'
1079 nrowa = n
1080 end if
1081 if (trans) then
1082 t = 'T'
1083 else
1084 t = 'N'
1085 end if
1086 if (present(err)) then
1087 errmgr => err
1088 else
1089 errmgr => deferr
1090 end if
1091
1092 ! Input Check
1093 flag = 0
1094 if (lside) then
1095 ! A is M-by-K, M >= K >= 0
1096 if (size(a, 1) /= m .or. size(a, 2) < k) flag = 3
1097 else
1098 ! A is N-by-K, N >= K >= 0
1099 if (size(a, 1) /= n .or. size(a, 2) < k) flag = 3
1100 end if
1101 if (flag /= 0) then
1102 ! ERROR: One of the input arrays is not sized correctly
1103 allocate(character(len = 256) :: errmsg)
1104 write(errmsg, 100) "Input number ", flag, &
1105 " is not sized correctly."
1106 call errmgr%report_error("mult_qr_mtx", trim(errmsg), &
1107 la_array_size_error)
1108 return
1109 end if
1110
1111 ! Workspace Query
1112 call dormqr(side, t, m, n, k, a, nrowa, tau, c, m, temp, -1, flag)
1113 lwork = int(temp(1), int32)
1114 if (present(olwork)) then
1115 olwork = lwork
1116 return
1117 end if
1118
1119 ! Local Memory Allocation
1120 if (present(work)) then
1121 if (size(work) < lwork) then
1122 ! ERROR: WORK not sized correctly
1123 call errmgr%report_error("mult_qr_mtx", &
1124 "Incorrectly sized input array WORK, argument 6.", &
1125 la_array_size_error)
1126 return
1127 end if
1128 wptr => work(1:lwork)
1129 else
1130 allocate(wrk(lwork), stat = istat)
1131 if (istat /= 0) then
1132 ! ERROR: Out of memory
1133 call errmgr%report_error("mult_qr_mtx", &
1134 "Insufficient memory available.", &
1135 la_out_of_memory_error)
1136 return
1137 end if
1138 wptr => wrk
1139 end if
1140
1141 ! Call DORMQR
1142 call dormqr(side, t, m, n, k, a, nrowa, tau, c, m, wptr, lwork, flag)
1143
1144 ! Formatting
1145100 format(a, i0, a)
1146 end subroutine
1147
1148! ------------------------------------------------------------------------------
1149 module subroutine mult_qr_mtx_cmplx(lside, trans, a, tau, c, work, olwork, err)
1150 ! Arguments
1151 logical, intent(in) :: lside, trans
1152 complex(real64), intent(in), dimension(:) :: tau
1153 complex(real64), intent(inout), dimension(:,:) :: a, c
1154 complex(real64), intent(out), target, dimension(:), optional :: work
1155 integer(int32), intent(out), optional :: olwork
1156 class(errors), intent(inout), optional, target :: err
1157
1158 ! Parameters
1159 complex(real64), parameter :: one = (1.0d0, 0.0d0)
1160
1161 ! Local Variables
1162 character :: side, t
1163 integer(int32) :: m, n, k, nrowa, istat, flag, lwork
1164 complex(real64), pointer, dimension(:) :: wptr
1165 complex(real64), allocatable, target, dimension(:) :: wrk
1166 complex(real64), dimension(1) :: temp
1167 class(errors), pointer :: errmgr
1168 type(errors), target :: deferr
1169 character(len = :), allocatable :: errmsg
1170
1171 ! Initialization
1172 m = size(c, 1)
1173 n = size(c, 2)
1174 k = size(tau)
1175 if (lside) then
1176 side = 'L'
1177 nrowa = m
1178 else
1179 side = 'R'
1180 nrowa = n
1181 end if
1182 if (trans) then
1183 t = 'C'
1184 else
1185 t = 'N'
1186 end if
1187 if (present(err)) then
1188 errmgr => err
1189 else
1190 errmgr => deferr
1191 end if
1192
1193 ! Input Check
1194 flag = 0
1195 if (lside) then
1196 ! A is M-by-K, M >= K >= 0
1197 if (size(a, 1) /= m .or. size(a, 2) < k) flag = 3
1198 else
1199 ! A is N-by-K, N >= K >= 0
1200 if (size(a, 1) /= n .or. size(a, 2) < k) flag = 3
1201 end if
1202 if (flag /= 0) then
1203 ! ERROR: One of the input arrays is not sized correctly
1204 allocate(character(len = 256) :: errmsg)
1205 write(errmsg, 100) "Input number ", flag, &
1206 " is not sized correctly."
1207 call errmgr%report_error("mult_qr_mtx_cmplx", trim(errmsg), &
1208 la_array_size_error)
1209 return
1210 end if
1211
1212 ! Workspace Query
1213 call zunmqr(side, t, m, n, k, a, nrowa, tau, c, m, temp, -1, flag)
1214 lwork = int(temp(1), int32)
1215 if (present(olwork)) then
1216 olwork = lwork
1217 return
1218 end if
1219
1220 ! Local Memory Allocation
1221 if (present(work)) then
1222 if (size(work) < lwork) then
1223 ! ERROR: WORK not sized correctly
1224 call errmgr%report_error("mult_qr_mtx_cmplx", &
1225 "Incorrectly sized input array WORK, argument 6.", &
1226 la_array_size_error)
1227 return
1228 end if
1229 wptr => work(1:lwork)
1230 else
1231 allocate(wrk(lwork), stat = istat)
1232 if (istat /= 0) then
1233 ! ERROR: Out of memory
1234 call errmgr%report_error("mult_qr_mtx_cmplx", &
1235 "Insufficient memory available.", &
1236 la_out_of_memory_error)
1237 return
1238 end if
1239 wptr => wrk
1240 end if
1241
1242 ! Call ZUNMQR
1243 call zunmqr(side, t, m, n, k, a, nrowa, tau, c, m, wptr, lwork, flag)
1244
1245 ! Formatting
1246100 format(a, i0, a)
1247 end subroutine
1248
1249! ------------------------------------------------------------------------------
1250 module subroutine mult_qr_vec(trans, a, tau, c, work, olwork, err)
1251 ! Arguments
1252 logical, intent(in) :: trans
1253 real(real64), intent(inout), dimension(:,:) :: a
1254 real(real64), intent(in), dimension(:) :: tau
1255 real(real64), intent(inout), dimension(:) :: c
1256 real(real64), intent(out), target, dimension(:), optional :: work
1257 integer(int32), intent(out), optional :: olwork
1258 class(errors), intent(inout), optional, target :: err
1259
1260 ! Parameters
1261 real(real64), parameter :: one = 1.0d0
1262
1263 ! Local Variables
1264 character :: side, t
1265 integer(int32) :: m, k, nrowa, istat, flag, lwork
1266 real(real64), pointer, dimension(:) :: wptr
1267 real(real64), allocatable, target, dimension(:) :: wrk
1268 real(real64), dimension(1) :: temp
1269 class(errors), pointer :: errmgr
1270 type(errors), target :: deferr
1271 character(len = :), allocatable :: errmsg
1272
1273 ! Initialization
1274 m = size(c)
1275 k = size(tau)
1276 side = 'L'
1277 nrowa = m
1278 if (trans) then
1279 t = 'T'
1280 else
1281 t = 'N'
1282 end if
1283 if (present(err)) then
1284 errmgr => err
1285 else
1286 errmgr => deferr
1287 end if
1288
1289 ! Input Check
1290 flag = 0
1291 if (size(a, 1) /= m .or. size(a, 2) < k) flag = 3
1292 if (flag /= 0) then
1293 ! ERROR: One of the input arrays is not sized correctly
1294 allocate(character(len = 256) :: errmsg)
1295 write(errmsg, 100) "Input number ", flag, &
1296 " is not sized correctly."
1297 call errmgr%report_error("mult_qr_vec", trim(errmsg), &
1298 la_array_size_error)
1299 return
1300 end if
1301
1302 ! Workspace Query
1303 call dormqr(side, t, m, 1, k, a, nrowa, tau, c, m, temp, -1, flag)
1304 lwork = int(temp(1), int32)
1305 if (present(olwork)) then
1306 olwork = lwork
1307 return
1308 end if
1309
1310 ! Local Memory Allocation
1311 if (present(work)) then
1312 if (size(work) < lwork) then
1313 ! ERROR: WORK not sized correctly
1314 call errmgr%report_error("mult_qr_vec", &
1315 "Incorrectly sized input array WORK, argument 6.", &
1316 la_array_size_error)
1317 return
1318 end if
1319 wptr => work(1:lwork)
1320 else
1321 allocate(wrk(lwork), stat = istat)
1322 if (istat /= 0) then
1323 ! ERROR: Out of memory
1324 call errmgr%report_error("mult_qr_vec", &
1325 "Insufficient memory available.", &
1326 la_out_of_memory_error)
1327 return
1328 end if
1329 wptr => wrk
1330 end if
1331
1332 ! Call DORMQR
1333 call dormqr(side, t, m, 1, k, a, nrowa, tau, c, m, wptr, lwork, flag)
1334
1335 ! Formatting
1336100 format(a, i0, a)
1337 end subroutine
1338
1339! ------------------------------------------------------------------------------
1340 module subroutine mult_qr_vec_cmplx(trans, a, tau, c, work, olwork, err)
1341 ! Arguments
1342 logical, intent(in) :: trans
1343 complex(real64), intent(inout), dimension(:,:) :: a
1344 complex(real64), intent(in), dimension(:) :: tau
1345 complex(real64), intent(inout), dimension(:) :: c
1346 complex(real64), intent(out), target, dimension(:), optional :: work
1347 integer(int32), intent(out), optional :: olwork
1348 class(errors), intent(inout), optional, target :: err
1349
1350 ! Parameters
1351 complex(real64), parameter :: one = (1.0d0, 0.0d0)
1352
1353 ! Local Variables
1354 character :: side, t
1355 integer(int32) :: m, k, nrowa, istat, flag, lwork
1356 complex(real64), pointer, dimension(:) :: wptr
1357 complex(real64), allocatable, target, dimension(:) :: wrk
1358 complex(real64), dimension(1) :: temp
1359 class(errors), pointer :: errmgr
1360 type(errors), target :: deferr
1361 character(len = :), allocatable :: errmsg
1362
1363 ! Initialization
1364 m = size(c)
1365 k = size(tau)
1366 side = 'L'
1367 nrowa = m
1368 if (trans) then
1369 t = 'C'
1370 else
1371 t = 'N'
1372 end if
1373 if (present(err)) then
1374 errmgr => err
1375 else
1376 errmgr => deferr
1377 end if
1378
1379 ! Input Check
1380 flag = 0
1381 if (size(a, 1) /= m .or. size(a, 2) < k) flag = 3
1382 if (flag /= 0) then
1383 ! ERROR: One of the input arrays is not sized correctly
1384 allocate(character(len = 256) :: errmsg)
1385 write(errmsg, 100) "Input number ", flag, &
1386 " is not sized correctly."
1387 call errmgr%report_error("mult_qr_vec", trim(errmsg), &
1388 la_array_size_error)
1389 return
1390 end if
1391
1392 ! Workspace Query
1393 call zunmqr(side, t, m, 1, k, a, nrowa, tau, c, m, temp, -1, flag)
1394 lwork = int(temp(1), int32)
1395 if (present(olwork)) then
1396 olwork = lwork
1397 return
1398 end if
1399
1400 ! Local Memory Allocation
1401 if (present(work)) then
1402 if (size(work) < lwork) then
1403 ! ERROR: WORK not sized correctly
1404 call errmgr%report_error("mult_qr_vec", &
1405 "Incorrectly sized input array WORK, argument 6.", &
1406 la_array_size_error)
1407 return
1408 end if
1409 wptr => work(1:lwork)
1410 else
1411 allocate(wrk(lwork), stat = istat)
1412 if (istat /= 0) then
1413 ! ERROR: Out of memory
1414 call errmgr%report_error("mult_qr_vec", &
1415 "Insufficient memory available.", &
1416 la_out_of_memory_error)
1417 return
1418 end if
1419 wptr => wrk
1420 end if
1421
1422 ! Call ZUNMQR
1423 call zunmqr(side, t, m, 1, k, a, nrowa, tau, c, m, wptr, lwork, flag)
1424
1425 ! Formatting
1426100 format(a, i0, a)
1427 end subroutine
1428
1429! ------------------------------------------------------------------------------
1430 module subroutine qr_rank1_update_dbl(q, r, u, v, work, err)
1431 ! Arguments
1432 real(real64), intent(inout), dimension(:,:) :: q, r
1433 real(real64), intent(inout), dimension(:) :: u, v
1434 real(real64), intent(out), target, optional, dimension(:) :: work
1435 class(errors), intent(inout), optional, target :: err
1436
1437 ! Local Variables
1438 logical :: full
1439 integer(int32) :: m, n, k, lwork, istat, flag
1440 real(real64), pointer, dimension(:) :: wptr
1441 real(real64), allocatable, target, dimension(:) :: wrk
1442 class(errors), pointer :: errmgr
1443 type(errors), target :: deferr
1444 character(len = :), allocatable :: errmsg
1445
1446 ! Initialization
1447 m = size(u, 1)
1448 n = size(r, 2)
1449 k = min(m, n)
1450 full = size(q, 2) == m
1451 lwork = 2 * k
1452 if (present(err)) then
1453 errmgr => err
1454 else
1455 errmgr => deferr
1456 end if
1457
1458 ! Input Check
1459 flag = 0
1460 if (m < n) then
1461 flag = 1
1462 else if (.not.full .and. size(q, 2) /= k) then
1463 flag = 1
1464 else if (size(r, 1) /= m) then
1465 flag = 2
1466 else if (size(u) /= m) then
1467 flag = 3
1468 else if (size(v) /= n) then
1469 flag = 4
1470 end if
1471 if (flag /= 0) then
1472 ! ERROR: One of the input arrays is not sized correctly
1473 allocate(character(len = 256) :: errmsg)
1474 write(errmsg, 100) "Input number ", flag, &
1475 " is not sized correctly."
1476 call errmgr%report_error("qr_rank1_update", trim(errmsg), &
1477 la_array_size_error)
1478 return
1479 end if
1480
1481 ! Local Memory Allocation
1482 if (present(work)) then
1483 if (size(work) < lwork) then
1484 ! ERROR: WORK not sized correctly
1485 call errmgr%report_error("qr_rank1_update", &
1486 "Incorrectly sized input array WORK, argument 5.", &
1487 la_array_size_error)
1488 return
1489 end if
1490 wptr => work(1:lwork)
1491 else
1492 allocate(wrk(lwork), stat = istat)
1493 if (istat /= 0) then
1494 ! ERROR: Out of memory
1495 call errmgr%report_error("qr_rank1_update", &
1496 "Insufficient memory available.", &
1497 la_out_of_memory_error)
1498 return
1499 end if
1500 wptr => wrk
1501 end if
1502
1503 ! Process
1504 call dqr1up(m, n, k, q, m, r, m, u, v, wptr)
1505
1506 ! Formatting
1507100 format(a, i0, a)
1508 end subroutine
1509
1510! ------------------------------------------------------------------------------
1511 module subroutine qr_rank1_update_cmplx(q, r, u, v, work, rwork, err)
1512 ! Arguments
1513 complex(real64), intent(inout), dimension(:,:) :: q, r
1514 complex(real64), intent(inout), dimension(:) :: u, v
1515 complex(real64), intent(out), target, optional, dimension(:) :: work
1516 real(real64), intent(out), target, optional, dimension(:) :: rwork
1517 class(errors), intent(inout), optional, target :: err
1518
1519 ! Local Variables
1520 logical :: full
1521 integer(int32) :: m, n, k, lwork, istat, flag, lrwork
1522 complex(real64), pointer, dimension(:) :: wptr
1523 complex(real64), allocatable, target, dimension(:) :: wrk
1524 real(real64), pointer, dimension(:) :: rwptr
1525 real(real64), allocatable, target, dimension(:) :: rwrk
1526 class(errors), pointer :: errmgr
1527 type(errors), target :: deferr
1528 character(len = :), allocatable :: errmsg
1529
1530 ! Initialization
1531 m = size(u, 1)
1532 n = size(r, 2)
1533 k = min(m, n)
1534 full = size(q, 2) == m
1535 lwork = k
1536 lrwork = k
1537 if (present(err)) then
1538 errmgr => err
1539 else
1540 errmgr => deferr
1541 end if
1542
1543 ! Input Check
1544 flag = 0
1545 if (m < n) then
1546 flag = 1
1547 else if (.not.full .and. size(q, 2) /= k) then
1548 flag = 1
1549 else if (size(r, 1) /= m) then
1550 flag = 2
1551 else if (size(u) /= m) then
1552 flag = 3
1553 else if (size(v) /= n) then
1554 flag = 4
1555 end if
1556 if (flag /= 0) then
1557 ! ERROR: One of the input arrays is not sized correctly
1558 allocate(character(len = 256) :: errmsg)
1559 write(errmsg, 100) "Input number ", flag, &
1560 " is not sized correctly."
1561 call errmgr%report_error("qr_rank1_update_cmplx", trim(errmsg), &
1562 la_array_size_error)
1563 return
1564 end if
1565
1566 ! Local Memory Allocation
1567 if (present(work)) then
1568 if (size(work) < lwork) then
1569 ! ERROR: WORK not sized correctly
1570 call errmgr%report_error("qr_rank1_update_cmplx", &
1571 "Incorrectly sized input array WORK, argument 5.", &
1572 la_array_size_error)
1573 return
1574 end if
1575 wptr => work(1:lwork)
1576 else
1577 allocate(wrk(lwork), stat = istat)
1578 if (istat /= 0) then
1579 ! ERROR: Out of memory
1580 call errmgr%report_error("qr_rank1_update_cmplx", &
1581 "Insufficient memory available.", &
1582 la_out_of_memory_error)
1583 return
1584 end if
1585 wptr => wrk
1586 end if
1587
1588 if (present(rwork)) then
1589 if (size(rwork) < lrwork) then
1590 ! ERROR: WORK not sized correctly
1591 call errmgr%report_error("qr_rank1_update_cmplx", &
1592 "Incorrectly sized input array RWORK, argument 6.", &
1593 la_array_size_error)
1594 return
1595 end if
1596 wptr => work(1:lrwork)
1597 else
1598 allocate(rwrk(lrwork), stat = istat)
1599 if (istat /= 0) then
1600 ! ERROR: Out of memory
1601 call errmgr%report_error("qr_rank1_update_cmplx", &
1602 "Insufficient memory available.", &
1603 la_out_of_memory_error)
1604 return
1605 end if
1606 rwptr => rwrk
1607 end if
1608
1609 ! Process
1610 call zqr1up(m, n, k, q, m, r, m, u, v, wptr, rwptr)
1611
1612 ! Formatting
1613100 format(a, i0, a)
1614 end subroutine
1615
1616! ******************************************************************************
1617! CHOLESKY FACTORIZATION
1618! ------------------------------------------------------------------------------
1619 module subroutine cholesky_factor_dbl(a, upper, err)
1620 ! Arguments
1621 real(real64), intent(inout), dimension(:,:) :: a
1622 logical, intent(in), optional :: upper
1623 class(errors), intent(inout), optional, target :: err
1624
1625 ! Parameters
1626 real(real64), parameter :: zero = 0.0d0
1627
1628 ! Local Variables
1629 character :: uplo
1630 integer(int32) :: i, n, flag
1631 class(errors), pointer :: errmgr
1632 type(errors), target :: deferr
1633 character(len = :), allocatable :: errmsg
1634
1635 ! Initialization
1636 n = size(a, 1)
1637 if (present(upper)) then
1638 if (upper) then
1639 uplo = 'U'
1640 else
1641 uplo = 'L'
1642 end if
1643 else
1644 uplo = 'U'
1645 end if
1646 if (present(err)) then
1647 errmgr => err
1648 else
1649 errmgr => deferr
1650 end if
1651
1652 ! Input Check
1653 if (size(a, 2) /= n) then
1654 ! ERROR: A must be square
1655 call errmgr%report_error("cholesky_factor", &
1656 "The input matrix must be square.", la_array_size_error)
1657 return
1658 end if
1659
1660 ! Process
1661 call dpotrf(uplo, n, a, n, flag)
1662 if (flag > 0) then
1663 ! ERROR: Matrix is not positive definite
1664 allocate(character(len = 256) :: errmsg)
1665 write(errmsg, 100) "The leading minor of order ", flag, &
1666 " is not positive definite."
1667 call errmgr%report_error("cholesky_factor", trim(errmsg), &
1668 la_matrix_format_error)
1669 end if
1670
1671 ! Zero out the non-used upper or lower diagonal
1672 if (uplo == 'U') then
1673 ! Zero out the lower
1674 do i = 1, n - 1
1675 a(i+1:n,i) = zero
1676 end do
1677 else
1678 ! Zero out the upper
1679 do i = 2, n
1680 a(1:i-1,i) = zero
1681 end do
1682 end if
1683
1684 ! Formatting
1685100 format(a, i0, a)
1686 end subroutine
1687
1688! ------------------------------------------------------------------------------
1689 module subroutine cholesky_factor_cmplx(a, upper, err)
1690 ! Arguments
1691 complex(real64), intent(inout), dimension(:,:) :: a
1692 logical, intent(in), optional :: upper
1693 class(errors), intent(inout), optional, target :: err
1694
1695 ! Parameters
1696 complex(real64), parameter :: zero = (0.0d0, 0.0d0)
1697
1698 ! Local Variables
1699 character :: uplo
1700 integer(int32) :: i, n, flag
1701 class(errors), pointer :: errmgr
1702 type(errors), target :: deferr
1703 character(len = :), allocatable :: errmsg
1704
1705 ! Initialization
1706 n = size(a, 1)
1707 if (present(upper)) then
1708 if (upper) then
1709 uplo = 'U'
1710 else
1711 uplo = 'L'
1712 end if
1713 else
1714 uplo = 'U'
1715 end if
1716 if (present(err)) then
1717 errmgr => err
1718 else
1719 errmgr => deferr
1720 end if
1721
1722 ! Input Check
1723 if (size(a, 2) /= n) then
1724 ! ERROR: A must be square
1725 call errmgr%report_error("cholesky_factor_cmplx", &
1726 "The input matrix must be square.", la_array_size_error)
1727 return
1728 end if
1729
1730 ! Process
1731 call zpotrf(uplo, n, a, n, flag)
1732 if (flag > 0) then
1733 ! ERROR: Matrix is not positive definite
1734 allocate(character(len = 256) :: errmsg)
1735 write(errmsg, 100) "The leading minor of order ", flag, &
1736 " is not positive definite."
1737 call errmgr%report_error("cholesky_factor_cmplx", trim(errmsg), &
1738 la_matrix_format_error)
1739 end if
1740
1741 ! Zero out the non-used upper or lower diagonal
1742 if (uplo == 'U') then
1743 ! Zero out the lower
1744 do i = 1, n - 1
1745 a(i+1:n,i) = zero
1746 end do
1747 else
1748 ! Zero out the upper
1749 do i = 2, n
1750 a(1:i-1,i) = zero
1751 end do
1752 end if
1753
1754 ! Formatting
1755100 format(a, i0, a)
1756 end subroutine
1757
1758! ------------------------------------------------------------------------------
1759 module subroutine cholesky_rank1_update_dbl(r, u, work, err)
1760 ! Arguments
1761 real(real64), intent(inout), dimension(:,:) :: r
1762 real(real64), intent(inout), dimension(:) :: u
1763 real(real64), intent(out), target, optional, dimension(:) :: work
1764 class(errors), intent(inout), optional, target :: err
1765
1766 ! Local Variables
1767 integer(int32) :: n, lwork, istat, flag
1768 real(real64), pointer, dimension(:) :: wptr
1769 real(real64), allocatable, target, dimension(:) :: wrk
1770 class(errors), pointer :: errmgr
1771 type(errors), target :: deferr
1772 character(len = :), allocatable :: errmsg
1773
1774 ! Initialization
1775 n = size(r, 1)
1776 lwork = n
1777 if (present(err)) then
1778 errmgr => err
1779 else
1780 errmgr => deferr
1781 end if
1782
1783 ! Input Check
1784 flag = 0
1785 if (size(r, 2) /= n) then
1786 flag = 1
1787 else if (size(u) /= n) then
1788 flag = 2
1789 end if
1790 if (flag /= 0) then
1791 allocate(character(len = 256) :: errmsg)
1792 write(errmsg, 100) "Input number ", flag, &
1793 " is not sized correctly."
1794 call errmgr%report_error("cholesky_rank1_update", trim(errmsg), &
1795 la_array_size_error)
1796 return
1797 end if
1798
1799 ! Local Memory Allocation
1800 if (present(work)) then
1801 if (size(work) < lwork) then
1802 ! ERROR: Workspace array is not sized correctly
1803 call errmgr%report_error("cholesky_rank1_update", &
1804 "The workspace array is too short.", &
1805 la_array_size_error)
1806 return
1807 end if
1808 wptr => work(1:lwork)
1809 else
1810 allocate(wrk(lwork), stat = istat)
1811 if (istat /= 0) then
1812 call errmgr%report_error("cholesky_rank1_update", &
1813 "Insufficient memory available.", &
1814 la_out_of_memory_error)
1815 return
1816 end if
1817 wptr => wrk
1818 end if
1819
1820 ! Process
1821 call dch1up(n, r, n, u, wptr)
1822
1823 ! Formatting
1824100 format(a, i0, a)
1825 end subroutine
1826
1827! ------------------------------------------------------------------------------
1828 module subroutine cholesky_rank1_update_cmplx(r, u, work, err)
1829 ! Arguments
1830 complex(real64), intent(inout), dimension(:,:) :: r
1831 complex(real64), intent(inout), dimension(:) :: u
1832 real(real64), intent(out), target, optional, dimension(:) :: work
1833 class(errors), intent(inout), optional, target :: err
1834
1835 ! Local Variables
1836 integer(int32) :: n, lwork, istat, flag
1837 real(real64), pointer, dimension(:) :: wptr
1838 real(real64), allocatable, target, dimension(:) :: wrk
1839 class(errors), pointer :: errmgr
1840 type(errors), target :: deferr
1841 character(len = :), allocatable :: errmsg
1842
1843 ! Initialization
1844 n = size(r, 1)
1845 lwork = n
1846 if (present(err)) then
1847 errmgr => err
1848 else
1849 errmgr => deferr
1850 end if
1851
1852 ! Input Check
1853 flag = 0
1854 if (size(r, 2) /= n) then
1855 flag = 1
1856 else if (size(u) /= n) then
1857 flag = 2
1858 end if
1859 if (flag /= 0) then
1860 allocate(character(len = 256) :: errmsg)
1861 write(errmsg, 100) "Input number ", flag, &
1862 " is not sized correctly."
1863 call errmgr%report_error("cholesky_rank1_update_cmplx", &
1864 trim(errmsg), &
1865 la_array_size_error)
1866 return
1867 end if
1868
1869 ! Local Memory Allocation
1870 if (present(work)) then
1871 if (size(work) < lwork) then
1872 ! ERROR: Workspace array is not sized correctly
1873 call errmgr%report_error("cholesky_rank1_update_cmplx", &
1874 "The workspace array is too short.", &
1875 la_array_size_error)
1876 return
1877 end if
1878 wptr => work(1:lwork)
1879 else
1880 allocate(wrk(lwork), stat = istat)
1881 if (istat /= 0) then
1882 call errmgr%report_error("cholesky_rank1_update", &
1883 "Insufficient memory available.", &
1884 la_out_of_memory_error)
1885 return
1886 end if
1887 wptr => wrk
1888 end if
1889
1890 ! Process
1891 call zch1up(n, r, n, u, wptr)
1892
1893 ! Formatting
1894100 format(a, i0, a)
1895 end subroutine
1896
1897! ------------------------------------------------------------------------------
1898 module subroutine cholesky_rank1_downdate_dbl(r, u, work, err)
1899 ! Arguments
1900 real(real64), intent(inout), dimension(:,:) :: r
1901 real(real64), intent(inout), dimension(:) :: u
1902 real(real64), intent(out), target, optional, dimension(:) :: work
1903 class(errors), intent(inout), optional, target :: err
1904
1905 ! Local Variables
1906 integer(int32) :: n, lwork, istat, flag
1907 real(real64), pointer, dimension(:) :: wptr
1908 real(real64), allocatable, target, dimension(:) :: wrk
1909 class(errors), pointer :: errmgr
1910 type(errors), target :: deferr
1911 character(len = :), allocatable :: errmsg
1912
1913 ! Initialization
1914 n = size(r, 1)
1915 lwork = n
1916 if (present(err)) then
1917 errmgr => err
1918 else
1919 errmgr => deferr
1920 end if
1921
1922 ! Input Check
1923 flag = 0
1924 if (size(r, 2) /= n) then
1925 flag = 1
1926 else if (size(u) /= n) then
1927 flag = 2
1928 end if
1929 if (flag /= 0) then
1930 allocate(character(len = 256) :: errmsg)
1931 write(errmsg, 100) "Input number ", flag, &
1932 " is not sized correctly."
1933 call errmgr%report_error("cholesky_rank1_downdate", trim(errmsg), &
1934 la_array_size_error)
1935 return
1936 end if
1937
1938 ! Local Memory Allocation
1939 if (present(work)) then
1940 if (size(work) < lwork) then
1941 ! ERROR: Workspace array is not sized correctly
1942 call errmgr%report_error("cholesky_rank1_downdate", &
1943 "The workspace array is too short.", &
1944 la_array_size_error)
1945 return
1946 end if
1947 wptr => work(1:lwork)
1948 else
1949 allocate(wrk(lwork), stat = istat)
1950 if (istat /= 0) then
1951 call errmgr%report_error("cholesky_rank1_downdate", &
1952 "Insufficient memory available.", &
1953 la_out_of_memory_error)
1954 return
1955 end if
1956 wptr => wrk
1957 end if
1958
1959 ! Process
1960 call dch1dn(n, r, n, u, wptr, flag)
1961 if (flag == 1) then
1962 ! ERROR: The matrix is not positive definite
1963 call errmgr%report_error("cholesky_rank1_downdate", &
1964 "The downdated matrix is not positive definite.", &
1965 la_matrix_format_error)
1966 else if (flag == 2) then
1967 ! ERROR: The matrix is singular
1968 call errmgr%report_error("cholesky_rank1_downdate", &
1969 "The input matrix is singular.", la_singular_matrix_error)
1970 end if
1971
1972 ! Formatting
1973100 format(a, i0, a)
1974 end subroutine
1975
1976! ------------------------------------------------------------------------------
1977 module subroutine cholesky_rank1_downdate_cmplx(r, u, work, err)
1978 ! Arguments
1979 complex(real64), intent(inout), dimension(:,:) :: r
1980 complex(real64), intent(inout), dimension(:) :: u
1981 real(real64), intent(out), target, optional, dimension(:) :: work
1982 class(errors), intent(inout), optional, target :: err
1983
1984 ! Local Variables
1985 integer(int32) :: n, lwork, istat, flag
1986 real(real64), pointer, dimension(:) :: wptr
1987 real(real64), allocatable, target, dimension(:) :: wrk
1988 class(errors), pointer :: errmgr
1989 type(errors), target :: deferr
1990 character(len = :), allocatable :: errmsg
1991
1992 ! Initialization
1993 n = size(r, 1)
1994 lwork = n
1995 if (present(err)) then
1996 errmgr => err
1997 else
1998 errmgr => deferr
1999 end if
2000
2001 ! Input Check
2002 flag = 0
2003 if (size(r, 2) /= n) then
2004 flag = 1
2005 else if (size(u) /= n) then
2006 flag = 2
2007 end if
2008 if (flag /= 0) then
2009 allocate(character(len = 256) :: errmsg)
2010 write(errmsg, 100) "Input number ", flag, &
2011 " is not sized correctly."
2012 call errmgr%report_error("cholesky_rank1_downdate_cmplx", &
2013 trim(errmsg), &
2014 la_array_size_error)
2015 return
2016 end if
2017
2018 ! Local Memory Allocation
2019 if (present(work)) then
2020 if (size(work) < lwork) then
2021 ! ERROR: Workspace array is not sized correctly
2022 call errmgr%report_error("cholesky_rank1_downdate_cmplx", &
2023 "The workspace array is too short.", &
2024 la_array_size_error)
2025 return
2026 end if
2027 wptr => work(1:lwork)
2028 else
2029 allocate(wrk(lwork), stat = istat)
2030 if (istat /= 0) then
2031 call errmgr%report_error("cholesky_rank1_downdate_cmplx", &
2032 "Insufficient memory available.", &
2033 la_out_of_memory_error)
2034 return
2035 end if
2036 wptr => wrk
2037 end if
2038
2039 ! Process
2040 call zch1dn(n, r, n, u, wptr, flag)
2041 if (flag == 1) then
2042 ! ERROR: The matrix is not positive definite
2043 call errmgr%report_error("cholesky_rank1_downdate_cmplx", &
2044 "The downdated matrix is not positive definite.", &
2045 la_matrix_format_error)
2046 else if (flag == 2) then
2047 ! ERROR: The matrix is singular
2048 call errmgr%report_error("cholesky_rank1_downdate_cmplx", &
2049 "The input matrix is singular.", la_singular_matrix_error)
2050 end if
2051
2052 ! Formatting
2053100 format(a, i0, a)
2054 end subroutine
2055
2056! ******************************************************************************
2057! RZ FACTORIZATION ROUTINES
2058! ------------------------------------------------------------------------------
2059 module subroutine rz_factor_dbl(a, tau, work, olwork, err)
2060 ! Arguments
2061 real(real64), intent(inout), dimension(:,:) :: a
2062 real(real64), intent(out), dimension(:) :: tau
2063 real(real64), intent(out), target, optional, dimension(:) :: work
2064 integer(int32), intent(out), optional :: olwork
2065 class(errors), intent(inout), optional, target :: err
2066
2067 ! Local Variables
2068 integer(int32) :: m, n, lwork, flag, istat
2069 real(real64), pointer, dimension(:) :: wptr
2070 real(real64), allocatable, target, dimension(:) :: wrk
2071 real(real64), dimension(1) :: temp
2072 class(errors), pointer :: errmgr
2073 type(errors), target :: deferr
2074 character(len = :), allocatable :: errmsg
2075
2076 ! Initialization
2077 m = size(a, 1)
2078 n = size(a, 2)
2079 if (present(err)) then
2080 errmgr => err
2081 else
2082 errmgr => deferr
2083 end if
2084
2085 ! Input Check
2086 flag = 0
2087 if (size(tau) /= m) then
2088 flag = 3
2089 end if
2090 if (flag /= 0) then
2091 ! ERROR: One of the input arrays is not sized correctly
2092 allocate(character(len = 256) :: errmsg)
2093 write(errmsg, 100) "Input number ", flag, &
2094 " is not sized correctly."
2095 call errmgr%report_error("rz_factor", trim(errmsg), &
2096 la_array_size_error)
2097 return
2098 end if
2099
2100 ! Workspace Query
2101 call dtzrzf(m, n, a, m, tau, temp, -1, flag)
2102 lwork = int(temp(1), int32)
2103 if (present(olwork)) then
2104 olwork = lwork
2105 return
2106 end if
2107
2108 ! Local Memory Allocation
2109 if (present(work)) then
2110 if (size(work) < lwork) then
2111 ! ERROR: WORK not sized correctly
2112 call errmgr%report_error("rz_factor", &
2113 "Incorrectly sized input array WORK, argument 3.", &
2114 la_array_size_error)
2115 return
2116 end if
2117 wptr => work(1:lwork)
2118 else
2119 allocate(wrk(lwork), stat = istat)
2120 if (istat /= 0) then
2121 ! ERROR: Out of memory
2122 call errmgr%report_error("rz_factor", &
2123 "Insufficient memory available.", &
2124 la_out_of_memory_error)
2125 return
2126 end if
2127 wptr => wrk
2128 end if
2129
2130 ! Call DTZRZF
2131 call dtzrzf(m, n, a, m, tau, wptr, lwork, flag)
2132
2133 ! Formatting
2134100 format(a, i0, a)
2135 end subroutine
2136
2137! ------------------------------------------------------------------------------
2138 module subroutine rz_factor_cmplx(a, tau, work, olwork, err)
2139 ! Arguments
2140 complex(real64), intent(inout), dimension(:,:) :: a
2141 complex(real64), intent(out), dimension(:) :: tau
2142 complex(real64), intent(out), target, optional, dimension(:) :: work
2143 integer(int32), intent(out), optional :: olwork
2144 class(errors), intent(inout), optional, target :: err
2145
2146 ! Local Variables
2147 integer(int32) :: m, n, lwork, flag, istat
2148 complex(real64), pointer, dimension(:) :: wptr
2149 complex(real64), allocatable, target, dimension(:) :: wrk
2150 complex(real64), dimension(1) :: temp
2151 class(errors), pointer :: errmgr
2152 type(errors), target :: deferr
2153 character(len = :), allocatable :: errmsg
2154
2155 ! Initialization
2156 m = size(a, 1)
2157 n = size(a, 2)
2158 if (present(err)) then
2159 errmgr => err
2160 else
2161 errmgr => deferr
2162 end if
2163
2164 ! Input Check
2165 flag = 0
2166 if (size(tau) /= m) then
2167 flag = 3
2168 end if
2169 if (flag /= 0) then
2170 ! ERROR: One of the input arrays is not sized correctly
2171 allocate(character(len = 256) :: errmsg)
2172 write(errmsg, 100) "Input number ", flag, &
2173 " is not sized correctly."
2174 call errmgr%report_error("rz_factor_cmplx", trim(errmsg), &
2175 la_array_size_error)
2176 return
2177 end if
2178
2179 ! Workspace Query
2180 call ztzrzf(m, n, a, m, tau, temp, -1, flag)
2181 lwork = int(temp(1), int32)
2182 if (present(olwork)) then
2183 olwork = lwork
2184 return
2185 end if
2186
2187 ! Local Memory Allocation
2188 if (present(work)) then
2189 if (size(work) < lwork) then
2190 ! ERROR: WORK not sized correctly
2191 call errmgr%report_error("rz_factor_cmplx", &
2192 "Incorrectly sized input array WORK, argument 3.", &
2193 la_array_size_error)
2194 return
2195 end if
2196 wptr => work(1:lwork)
2197 else
2198 allocate(wrk(lwork), stat = istat)
2199 if (istat /= 0) then
2200 ! ERROR: Out of memory
2201 call errmgr%report_error("rz_factor_cmplx", &
2202 "Insufficient memory available.", &
2203 la_out_of_memory_error)
2204 return
2205 end if
2206 wptr => wrk
2207 end if
2208
2209 ! Call ZTZRZF
2210 call ztzrzf(m, n, a, m, tau, wptr, lwork, flag)
2211
2212 ! Formatting
2213100 format(a, i0, a)
2214 end subroutine
2215
2216! ------------------------------------------------------------------------------
2217 module subroutine mult_rz_mtx(lside, trans, l, a, tau, c, work, olwork, err)
2218 ! Arguments
2219 logical, intent(in) :: lside, trans
2220 integer(int32), intent(in) :: l
2221 real(real64), intent(inout), dimension(:,:) :: a, c
2222 real(real64), intent(in), dimension(:) :: tau
2223 real(real64), intent(out), target, optional, dimension(:) :: work
2224 integer(int32), intent(out), optional :: olwork
2225 class(errors), intent(inout), optional, target :: err
2226
2227 ! Local Variables
2228 character :: side, t
2229 integer(int32) :: m, n, k, lwork, flag, istat, lda
2230 real(real64), pointer, dimension(:) :: wptr
2231 real(real64), allocatable, target, dimension(:) :: wrk
2232 real(real64), dimension(1) :: temp
2233 class(errors), pointer :: errmgr
2234 type(errors), target :: deferr
2235 character(len = :), allocatable :: errmsg
2236
2237 ! Initialization
2238 m = size(c, 1)
2239 n = size(c, 2)
2240 k = size(tau)
2241 lda = size(a, 1)
2242 if (lside) then
2243 side = 'L'
2244 else
2245 side = 'R'
2246 end if
2247 if (trans) then
2248 t = 'T'
2249 else
2250 t = 'N'
2251 end if
2252 if (present(err)) then
2253 errmgr => err
2254 else
2255 errmgr => deferr
2256 end if
2257
2258 ! Input Check
2259 flag = 0
2260 if (lside) then
2261 if (l > m .or. l < 0) then
2262 flag = 3
2263 else if (k > m) then
2264 flag = 5
2265 else if (size(a, 1) < k .or. size(a, 2) /= m) then
2266 flag = 4
2267 end if
2268 else
2269 if (l > n .or. l < 0) then
2270 flag = 3
2271 else if (k > n) then
2272 flag = 5
2273 else if (size(a, 1) < k .or. size(a, 2) /= n) then
2274 flag = 4
2275 end if
2276 end if
2277 if (flag /= 0) then
2278 ! ERROR: One of the input arrays is not sized correctly
2279 allocate(character(len = 256) :: errmsg)
2280 write(errmsg, 100) "Input number ", flag, &
2281 " is not sized correctly."
2282 call errmgr%report_error("mult_rz_mtx", trim(errmsg), &
2283 la_array_size_error)
2284 return
2285 end if
2286
2287 ! Workspace Query
2288 call dormrz(side, t, m, n, k, l, a, lda, tau, c, m, temp, -1, flag)
2289 lwork = int(temp(1), int32)
2290 if (present(olwork)) then
2291 olwork = lwork
2292 return
2293 end if
2294
2295 ! Local Memory Allocation
2296 if (present(work)) then
2297 if (size(work) < lwork) then
2298 ! ERROR: WORK not sized correctly
2299 call errmgr%report_error("mult_rz_mtx", &
2300 "Incorrectly sized input array WORK, argument 7.", &
2301 la_array_size_error)
2302 return
2303 end if
2304 wptr => work(1:lwork)
2305 else
2306 allocate(wrk(lwork), stat = istat)
2307 if (istat /= 0) then
2308 ! ERROR: Out of memory
2309 call errmgr%report_error("mult_rz_mtx", &
2310 "Insufficient memory available.", &
2311 la_out_of_memory_error)
2312 return
2313 end if
2314 wptr => wrk
2315 end if
2316
2317 ! Call DORMRZ
2318 call dormrz(side, t, m, n, k, l, a, lda, tau, c, m, wptr, lwork, flag)
2319
2320 ! Formatting
2321100 format(a, i0, a)
2322 end subroutine
2323
2324! ------------------------------------------------------------------------------
2325 module subroutine mult_rz_mtx_cmplx(lside, trans, l, a, tau, c, work, olwork, err)
2326 ! Arguments
2327 logical, intent(in) :: lside, trans
2328 integer(int32), intent(in) :: l
2329 complex(real64), intent(inout), dimension(:,:) :: a, c
2330 complex(real64), intent(in), dimension(:) :: tau
2331 complex(real64), intent(out), target, optional, dimension(:) :: work
2332 integer(int32), intent(out), optional :: olwork
2333 class(errors), intent(inout), optional, target :: err
2334
2335 ! Local Variables
2336 character :: side, t
2337 integer(int32) :: m, n, k, lwork, flag, istat, lda
2338 complex(real64), pointer, dimension(:) :: wptr
2339 complex(real64), allocatable, target, dimension(:) :: wrk
2340 complex(real64), dimension(1) :: temp
2341 class(errors), pointer :: errmgr
2342 type(errors), target :: deferr
2343 character(len = :), allocatable :: errmsg
2344
2345 ! Initialization
2346 m = size(c, 1)
2347 n = size(c, 2)
2348 k = size(tau)
2349 lda = size(a, 1)
2350 if (lside) then
2351 side = 'L'
2352 else
2353 side = 'R'
2354 end if
2355 if (trans) then
2356 t = 'C'
2357 else
2358 t = 'N'
2359 end if
2360 if (present(err)) then
2361 errmgr => err
2362 else
2363 errmgr => deferr
2364 end if
2365
2366 ! Input Check
2367 flag = 0
2368 if (lside) then
2369 if (l > m .or. l < 0) then
2370 flag = 3
2371 else if (k > m) then
2372 flag = 5
2373 else if (size(a, 1) < k .or. size(a, 2) /= m) then
2374 flag = 4
2375 end if
2376 else
2377 if (l > n .or. l < 0) then
2378 flag = 3
2379 else if (k > n) then
2380 flag = 5
2381 else if (size(a, 1) < k .or. size(a, 2) /= n) then
2382 flag = 4
2383 end if
2384 end if
2385 if (flag /= 0) then
2386 ! ERROR: One of the input arrays is not sized correctly
2387 allocate(character(len = 256) :: errmsg)
2388 write(errmsg, 100) "Input number ", flag, &
2389 " is not sized correctly."
2390 call errmgr%report_error("mult_rz_mtx_cmplx", trim(errmsg), &
2391 la_array_size_error)
2392 return
2393 end if
2394
2395 ! Workspace Query
2396 call zunmrz(side, t, m, n, k, l, a, lda, tau, c, m, temp, -1, flag)
2397 lwork = int(temp(1), int32)
2398 if (present(olwork)) then
2399 olwork = lwork
2400 return
2401 end if
2402
2403 ! Local Memory Allocation
2404 if (present(work)) then
2405 if (size(work) < lwork) then
2406 ! ERROR: WORK not sized correctly
2407 call errmgr%report_error("mult_rz_mtx_cmplx", &
2408 "Incorrectly sized input array WORK, argument 7.", &
2409 la_array_size_error)
2410 return
2411 end if
2412 wptr => work(1:lwork)
2413 else
2414 allocate(wrk(lwork), stat = istat)
2415 if (istat /= 0) then
2416 ! ERROR: Out of memory
2417 call errmgr%report_error("mult_rz_mtx_cmplx", &
2418 "Insufficient memory available.", &
2419 la_out_of_memory_error)
2420 return
2421 end if
2422 wptr => wrk
2423 end if
2424
2425 ! Call ZUNMRZ
2426 call zunmrz(side, t, m, n, k, l, a, lda, tau, c, m, wptr, lwork, flag)
2427
2428 ! Formatting
2429100 format(a, i0, a)
2430 end subroutine
2431
2432! ------------------------------------------------------------------------------
2433 module subroutine mult_rz_vec(trans, l, a, tau, c, work, olwork, err)
2434 ! Arguments
2435 logical, intent(in) :: trans
2436 integer(int32), intent(in) :: l
2437 real(real64), intent(inout), dimension(:,:) :: a
2438 real(real64), intent(in), dimension(:) :: tau
2439 real(real64), intent(inout), dimension(:) :: c
2440 real(real64), intent(out), target, optional, dimension(:) :: work
2441 integer(int32), intent(out), optional :: olwork
2442 class(errors), intent(inout), optional, target :: err
2443
2444 ! Local Variables
2445 character :: side, t
2446 integer(int32) :: m, k, lwork, flag, istat, lda
2447 real(real64), pointer, dimension(:) :: wptr
2448 real(real64), allocatable, target, dimension(:) :: wrk
2449 real(real64), dimension(1) :: temp
2450 class(errors), pointer :: errmgr
2451 type(errors), target :: deferr
2452 character(len = :), allocatable :: errmsg
2453
2454 ! Initialization
2455 m = size(c)
2456 k = size(tau)
2457 lda = size(a, 1)
2458 side = 'L'
2459 if (trans) then
2460 t = 'T'
2461 else
2462 t = 'N'
2463 end if
2464 if (present(err)) then
2465 errmgr => err
2466 else
2467 errmgr => deferr
2468 end if
2469
2470 ! Input Check
2471 flag = 0
2472 if (l > m .or. l < 0) then
2473 flag = 2
2474 else if (k > m) then
2475 flag = 4
2476 else if (size(a, 1) < k .or. size(a, 2) /= m) then
2477 flag = 3
2478 end if
2479 if (flag /= 0) then
2480 ! ERROR: One of the input arrays is not sized correctly
2481 allocate(character(len = 256) :: errmsg)
2482 write(errmsg, 100) "Input number ", flag, &
2483 " is not sized correctly."
2484 call errmgr%report_error("mult_rz_vec", trim(errmsg), &
2485 la_array_size_error)
2486 return
2487 end if
2488
2489 ! Workspace Query
2490 call dormrz(side, t, m, 1, k, l, a, lda, tau, c, m, temp, -1, flag)
2491 lwork = int(temp(1), int32)
2492 if (present(olwork)) then
2493 olwork = lwork
2494 return
2495 end if
2496
2497 ! Local Memory Allocation
2498 if (present(work)) then
2499 if (size(work) < lwork) then
2500 ! ERROR: WORK not sized correctly
2501 call errmgr%report_error("mult_rz_vec", &
2502 "Incorrectly sized input array WORK, argument 6.", &
2503 la_array_size_error)
2504 return
2505 end if
2506 wptr => work(1:lwork)
2507 else
2508 allocate(wrk(lwork), stat = istat)
2509 if (istat /= 0) then
2510 ! ERROR: Out of memory
2511 call errmgr%report_error("mult_rz_vec", &
2512 "Insufficient memory available.", &
2513 la_out_of_memory_error)
2514 return
2515 end if
2516 wptr => wrk
2517 end if
2518
2519 ! Call DORMRZ
2520 call dormrz(side, t, m, 1, k, l, a, lda, tau, c, m, wptr, lwork, flag)
2521
2522 ! Formatting
2523100 format(a, i0, a)
2524 end subroutine
2525
2526! ------------------------------------------------------------------------------
2527 module subroutine mult_rz_vec_cmplx(trans, l, a, tau, c, work, olwork, err)
2528 ! Arguments
2529 logical, intent(in) :: trans
2530 integer(int32), intent(in) :: l
2531 complex(real64), intent(inout), dimension(:,:) :: a
2532 complex(real64), intent(in), dimension(:) :: tau
2533 complex(real64), intent(inout), dimension(:) :: c
2534 complex(real64), intent(out), target, optional, dimension(:) :: work
2535 integer(int32), intent(out), optional :: olwork
2536 class(errors), intent(inout), optional, target :: err
2537
2538 ! Local Variables
2539 character :: side, t
2540 integer(int32) :: m, k, lwork, flag, istat, lda
2541 complex(real64), pointer, dimension(:) :: wptr
2542 complex(real64), allocatable, target, dimension(:) :: wrk
2543 complex(real64), dimension(1) :: temp
2544 class(errors), pointer :: errmgr
2545 type(errors), target :: deferr
2546 character(len = :), allocatable :: errmsg
2547
2548 ! Initialization
2549 m = size(c)
2550 k = size(tau)
2551 lda = size(a, 1)
2552 side = 'L'
2553 if (trans) then
2554 t = 'C'
2555 else
2556 t = 'N'
2557 end if
2558 if (present(err)) then
2559 errmgr => err
2560 else
2561 errmgr => deferr
2562 end if
2563
2564 ! Input Check
2565 flag = 0
2566 if (l > m .or. l < 0) then
2567 flag = 2
2568 else if (k > m) then
2569 flag = 4
2570 else if (size(a, 1) < k .or. size(a, 2) /= m) then
2571 flag = 3
2572 end if
2573 if (flag /= 0) then
2574 ! ERROR: One of the input arrays is not sized correctly
2575 allocate(character(len = 256) :: errmsg)
2576 write(errmsg, 100) "Input number ", flag, &
2577 " is not sized correctly."
2578 call errmgr%report_error("mult_rz_vec_cmplx", trim(errmsg), &
2579 la_array_size_error)
2580 return
2581 end if
2582
2583 ! Workspace Query
2584 call zunmrz(side, t, m, 1, k, l, a, lda, tau, c, m, temp, -1, flag)
2585 lwork = int(temp(1), int32)
2586 if (present(olwork)) then
2587 olwork = lwork
2588 return
2589 end if
2590
2591 ! Local Memory Allocation
2592 if (present(work)) then
2593 if (size(work) < lwork) then
2594 ! ERROR: WORK not sized correctly
2595 call errmgr%report_error("mult_rz_vec_cmplx", &
2596 "Incorrectly sized input array WORK, argument 6.", &
2597 la_array_size_error)
2598 return
2599 end if
2600 wptr => work(1:lwork)
2601 else
2602 allocate(wrk(lwork), stat = istat)
2603 if (istat /= 0) then
2604 ! ERROR: Out of memory
2605 call errmgr%report_error("mult_rz_vec_cmplx", &
2606 "Insufficient memory available.", &
2607 la_out_of_memory_error)
2608 return
2609 end if
2610 wptr => wrk
2611 end if
2612
2613 ! Call ZUNMRZ
2614 call zunmrz(side, t, m, 1, k, l, a, lda, tau, c, m, wptr, lwork, flag)
2615
2616 ! Formatting
2617100 format(a, i0, a)
2618 end subroutine
2619
2620! ******************************************************************************
2621! SVD ROUTINES
2622! ------------------------------------------------------------------------------
2623 module subroutine svd_dbl(a, s, u, vt, work, olwork, err)
2624 ! Arguments
2625 real(real64), intent(inout), dimension(:,:) :: a
2626 real(real64), intent(out), dimension(:) :: s
2627 real(real64), intent(out), optional, dimension(:,:) :: u, vt
2628 real(real64), intent(out), target, optional, dimension(:) :: work
2629 integer(int32), intent(out), optional :: olwork
2630 class(errors), intent(inout), optional, target :: err
2631
2632 ! Local Variables
2633 character :: jobu, jobvt
2634 integer(int32) :: m, n, mn, istat, lwork, flag
2635 real(real64), pointer, dimension(:) :: wptr
2636 real(real64), allocatable, target, dimension(:) :: wrk
2637 real(real64), dimension(1) :: temp
2638 class(errors), pointer :: errmgr
2639 type(errors), target :: deferr
2640 character(len = :), allocatable :: errmsg
2641
2642 ! Initialization
2643 m = size(a, 1)
2644 n = size(a, 2)
2645 mn = min(m, n)
2646 if (present(u)) then
2647 if (size(u, 2) == m) then
2648 jobu = 'A'
2649 else if (size(u, 2) == mn) then
2650 jobu = 'S'
2651 end if
2652 else
2653 jobu = 'N'
2654 end if
2655 if (present(vt)) then
2656 jobvt = 'A'
2657 else
2658 jobvt = 'N'
2659 end if
2660 if (present(err)) then
2661 errmgr => err
2662 else
2663 errmgr => deferr
2664 end if
2665
2666 ! Input Check
2667 flag = 0
2668 if (size(s) /= mn) then
2669 flag = 2
2670 else if (present(u)) then
2671 if (size(u, 1) /= m) flag = 3
2672 if (size(u, 2) /= m .and. size(u, 2) /= mn) flag = 3
2673 else if (present(vt)) then
2674 if (size(vt, 1) /= n .or. size(vt, 2) /= n) flag = 4
2675 end if
2676 if (flag /= 0) then
2677 ! ERROR: One of the input arrays is not sized correctly
2678 allocate(character(len = 256) :: errmsg)
2679 write(errmsg, 100) "Input number ", flag, &
2680 " is not sized correctly."
2681 call errmgr%report_error("svd", trim(errmsg), &
2682 la_array_size_error)
2683 return
2684 end if
2685
2686 ! Workspace Query
2687 call dgesvd(jobu, jobvt, m, n, a, m, s, temp, m, temp, n, temp, -1, &
2688 flag)
2689 lwork = int(temp(1), int32)
2690 if (present(olwork)) then
2691 olwork = lwork
2692 return
2693 end if
2694
2695 ! Local Memory Allocation
2696 if (present(work)) then
2697 if (size(work) < lwork) then
2698 ! ERROR: WORK not sized correctly
2699 call errmgr%report_error("svd", &
2700 "Incorrectly sized input array WORK, argument 5.", &
2701 la_array_size_error)
2702 return
2703 end if
2704 wptr => work(1:lwork)
2705 else
2706 allocate(wrk(lwork), stat = istat)
2707 if (istat /= 0) then
2708 ! ERROR: Out of memory
2709 call errmgr%report_error("svd", &
2710 "Insufficient memory available.", &
2711 la_out_of_memory_error)
2712 return
2713 end if
2714 wptr => wrk
2715 end if
2716
2717 ! Call DGESVD
2718 if (present(u) .and. present(vt)) then
2719 call dgesvd(jobu, jobvt, m, n, a, m, s, u, m, vt, n, wptr, lwork, &
2720 flag)
2721 else if (present(u) .and. .not.present(vt)) then
2722 call dgesvd(jobu, jobvt, m, n, a, m, s, u, m, temp, n, wptr, &
2723 lwork, flag)
2724 else if (.not.present(u) .and. present(vt)) then
2725 call dgesvd(jobu, jobvt, m, n, a, m, s, temp, m, vt, n, wptr, &
2726 lwork, flag)
2727 else
2728 call dgesvd(jobu, jobvt, m, n, a, m, s, temp, m, temp, n, wptr, &
2729 lwork, flag)
2730 end if
2731
2732 ! Check for convergence
2733 if (flag > 0) then
2734 allocate(character(len = 256) :: errmsg)
2735 write(errmsg, 101) flag, " superdiagonals could not " // &
2736 "converge to zero as part of the QR iteration process."
2737 call errmgr%report_warning("svd", errmsg, la_convergence_error)
2738 end if
2739
2740 ! Formatting
2741100 format(a, i0, a)
2742101 format(i0, a)
2743 end subroutine
2744
2745! ------------------------------------------------------------------------------
2746 module subroutine svd_cmplx(a, s, u, vt, work, olwork, rwork, err)
2747 ! Arguments
2748 complex(real64), intent(inout), dimension(:,:) :: a
2749 real(real64), intent(out), dimension(:) :: s
2750 complex(real64), intent(out), optional, dimension(:,:) :: u, vt
2751 complex(real64), intent(out), target, optional, dimension(:) :: work
2752 integer(int32), intent(out), optional :: olwork
2753 real(real64), intent(out), target, optional, dimension(:) :: rwork
2754 class(errors), intent(inout), optional, target :: err
2755
2756 ! Local Variables
2757 character :: jobu, jobvt
2758 integer(int32) :: m, n, mn, istat, lwork, flag, lrwork
2759 complex(real64), pointer, dimension(:) :: wptr
2760 complex(real64), allocatable, target, dimension(:) :: wrk
2761 complex(real64), dimension(1) :: temp
2762 real(real64), dimension(1) :: rtemp
2763 real(real64), pointer, dimension(:) :: rwptr
2764 real(real64), allocatable, target, dimension(:) :: rwrk
2765 class(errors), pointer :: errmgr
2766 type(errors), target :: deferr
2767 character(len = :), allocatable :: errmsg
2768
2769 ! Initialization
2770 m = size(a, 1)
2771 n = size(a, 2)
2772 mn = min(m, n)
2773 lrwork = 5 * mn
2774 if (present(u)) then
2775 if (size(u, 2) == m) then
2776 jobu = 'A'
2777 else if (size(u, 2) == mn) then
2778 jobu = 'S'
2779 end if
2780 else
2781 jobu = 'N'
2782 end if
2783 if (present(vt)) then
2784 jobvt = 'A'
2785 else
2786 jobvt = 'N'
2787 end if
2788 if (present(err)) then
2789 errmgr => err
2790 else
2791 errmgr => deferr
2792 end if
2793
2794 ! Input Check
2795 flag = 0
2796 if (size(s) /= mn) then
2797 flag = 2
2798 else if (present(u)) then
2799 if (size(u, 1) /= m) flag = 3
2800 if (size(u, 2) /= m .and. size(u, 2) /= mn) flag = 3
2801 else if (present(vt)) then
2802 if (size(vt, 1) /= n .or. size(vt, 2) /= n) flag = 4
2803 end if
2804 if (flag /= 0) then
2805 ! ERROR: One of the input arrays is not sized correctly
2806 allocate(character(len = 256) :: errmsg)
2807 write(errmsg, 100) "Input number ", flag, &
2808 " is not sized correctly."
2809 call errmgr%report_error("svd_cmplx", trim(errmsg), &
2810 la_array_size_error)
2811 return
2812 end if
2813
2814 ! Workspace Query
2815 call zgesvd(jobu, jobvt, m, n, a, m, s, temp, m, temp, n, temp, -1, &
2816 rtemp, flag)
2817 lwork = int(temp(1), int32)
2818 if (present(olwork)) then
2819 olwork = lwork
2820 return
2821 end if
2822
2823 ! Local Memory Allocation
2824 if (present(work)) then
2825 if (size(work) < lwork) then
2826 ! ERROR: WORK not sized correctly
2827 call errmgr%report_error("svd_cmplx", &
2828 "Incorrectly sized input array WORK, argument 5.", &
2829 la_array_size_error)
2830 return
2831 end if
2832 wptr => work(1:lwork)
2833 else
2834 allocate(wrk(lwork), stat = istat)
2835 if (istat /= 0) then
2836 ! ERROR: Out of memory
2837 call errmgr%report_error("svd_cmplx", &
2838 "Insufficient memory available.", &
2839 la_out_of_memory_error)
2840 return
2841 end if
2842 wptr => wrk
2843 end if
2844
2845 if (present(rwork)) then
2846 if (size(rwork) < lrwork) then
2847 ! ERROR: RWORK not sized correctly
2848 call errmgr%report_error("svd_cmplx", &
2849 "Incorrectly sized input array RWORK, argument 7.", &
2850 la_array_size_error)
2851 end if
2852 rwptr => rwork(1:lrwork)
2853 else
2854 allocate(rwrk(lrwork), stat = istat)
2855 if (istat /= 0) then
2856 ! ERROR: Out of memory
2857 call errmgr%report_error("svd_cmplx", &
2858 "Insufficient memory available.", &
2859 la_out_of_memory_error)
2860 return
2861 end if
2862 rwptr => rwrk
2863 end if
2864
2865 ! Call ZGESVD
2866 if (present(u) .and. present(vt)) then
2867 call zgesvd(jobu, jobvt, m, n, a, m, s, u, m, vt, n, wptr, lwork, &
2868 rwptr, flag)
2869 else if (present(u) .and. .not.present(vt)) then
2870 call zgesvd(jobu, jobvt, m, n, a, m, s, u, m, temp, n, wptr, &
2871 lwork, rwptr, flag)
2872 else if (.not.present(u) .and. present(vt)) then
2873 call zgesvd(jobu, jobvt, m, n, a, m, s, temp, m, vt, n, wptr, &
2874 lwork, rwptr, flag)
2875 else
2876 call zgesvd(jobu, jobvt, m, n, a, m, s, temp, m, temp, n, wptr, &
2877 lwork, rwptr, flag)
2878 end if
2879
2880 ! Check for convergence
2881 if (flag > 0) then
2882 allocate(character(len = 256) :: errmsg)
2883 write(errmsg, 101) flag, " superdiagonals could not " // &
2884 "converge to zero as part of the QR iteration process."
2885 call errmgr%report_warning("svd_cmplx", errmsg, &
2886 la_convergence_error)
2887 end if
2888
2889 ! Formatting
2890100 format(a, i0, a)
2891101 format(i0, a)
2892 end subroutine
2893
2894! ******************************************************************************
2895! LQ FACTORIZATION
2896! ------------------------------------------------------------------------------
2897 module subroutine lq_factor_no_pivot(a, tau, work, olwork, err)
2898 ! Arguments
2899 real(real64), intent(inout), dimension(:,:) :: a
2900 real(real64), intent(out), dimension(:) :: tau
2901 real(real64), intent(out), target, dimension(:), optional :: work
2902 integer(int32), intent(out), optional :: olwork
2903 class(errors), intent(inout), optional, target :: err
2904
2905 ! Local Variables
2906 integer(int32) :: m, n, mn, istat, lwork, flag
2907 real(real64), dimension(1) :: temp
2908 real(real64), pointer, dimension(:) :: wptr
2909 real(real64), allocatable, target, dimension(:) :: wrk
2910 class(errors), pointer :: errmgr
2911 type(errors), target :: deferr
2912
2913 ! Initialization
2914 m = size(a, 1)
2915 n = size(a, 2)
2916 mn = min(m, n)
2917 if (present(err)) then
2918 errmgr => err
2919 else
2920 errmgr => deferr
2921 end if
2922
2923 ! Input Check
2924 if (size(tau) /= mn) then
2925 ! ERROR: TAU not sized correctly
2926 call errmgr%report_error("lq_factor_no_pivot", &
2927 "Incorrectly sized input array TAU, argument 2.", &
2928 la_array_size_error)
2929 return
2930 end if
2931
2932 ! Workspace Query
2933 call dgelqf(m, n, a, m, tau, temp, -1, flag)
2934 lwork = int(temp(1), int32)
2935 if (present(olwork)) then
2936 olwork = lwork
2937 return
2938 end if
2939
2940 ! Local Memory Allocation
2941 if (present(work)) then
2942 if (size(work) < lwork) then
2943 ! ERROR: WORK not sized correctly
2944 call errmgr%report_error("lq_factor_no_pivot", &
2945 "Incorrectly sized input array WORK, argument 3.", &
2946 la_array_size_error)
2947 return
2948 end if
2949 wptr => work(1:lwork)
2950 else
2951 allocate(wrk(lwork), stat = istat)
2952 if (istat /= 0) then
2953 ! ERROR: Out of memory
2954 call errmgr%report_error("lq_factor_no_pivot", &
2955 "Insufficient memory available.", &
2956 la_out_of_memory_error)
2957 return
2958 end if
2959 wptr => wrk
2960 end if
2961
2962 ! Call DGELQF
2963 call dgelqf(m, n, a, m, tau, wptr, lwork, flag)
2964 end subroutine
2965
2966! ------------------------------------------------------------------------------
2967 module subroutine lq_factor_no_pivot_cmplx(a, tau, work, olwork, err)
2968 ! Arguments
2969 complex(real64), intent(inout), dimension(:,:) :: a
2970 complex(real64), intent(out), dimension(:) :: tau
2971 complex(real64), intent(out), target, dimension(:), optional :: work
2972 integer(int32), intent(out), optional :: olwork
2973 class(errors), intent(inout), optional, target :: err
2974
2975 ! Local Variables
2976 integer(int32) :: m, n, mn, istat, lwork, flag
2977 complex(real64), dimension(1) :: temp
2978 complex(real64), pointer, dimension(:) :: wptr
2979 complex(real64), allocatable, target, dimension(:) :: wrk
2980 class(errors), pointer :: errmgr
2981 type(errors), target :: deferr
2982
2983 ! Initialization
2984 m = size(a, 1)
2985 n = size(a, 2)
2986 mn = min(m, n)
2987 if (present(err)) then
2988 errmgr => err
2989 else
2990 errmgr => deferr
2991 end if
2992
2993 ! Input Check
2994 if (size(tau) /= mn) then
2995 ! ERROR: TAU not sized correctly
2996 call errmgr%report_error("lq_factor_no_pivot_cmplx", &
2997 "Incorrectly sized input array TAU, argument 2.", &
2998 la_array_size_error)
2999 return
3000 end if
3001
3002 ! Workspace Query
3003 call zgelqf(m, n, a, m, tau, temp, -1, flag)
3004 lwork = int(temp(1), int32)
3005 if (present(olwork)) then
3006 olwork = lwork
3007 return
3008 end if
3009
3010 ! Local Memory Allocation
3011 if (present(work)) then
3012 if (size(work) < lwork) then
3013 ! ERROR: WORK not sized correctly
3014 call errmgr%report_error("lq_factor_no_pivot_cmplx", &
3015 "Incorrectly sized input array WORK, argument 3.", &
3016 la_array_size_error)
3017 return
3018 end if
3019 wptr => work(1:lwork)
3020 else
3021 allocate(wrk(lwork), stat = istat)
3022 if (istat /= 0) then
3023 ! ERROR: Out of memory
3024 call errmgr%report_error("lq_factor_no_pivot_cmplx", &
3025 "Insufficient memory available.", &
3026 la_out_of_memory_error)
3027 return
3028 end if
3029 wptr => wrk
3030 end if
3031
3032 ! Call ZGELQF
3033 call zgelqf(m, n, a, m, tau, wptr, lwork, flag)
3034 end subroutine
3035
3036! ------------------------------------------------------------------------------
3037 module subroutine form_lq_no_pivot(l, tau, q, work, olwork, err)
3038 ! Arguments
3039 real(real64), intent(inout), dimension(:,:) :: l
3040 real(real64), intent(in), dimension(:) :: tau
3041 real(real64), intent(out), dimension(:,:) :: q
3042 real(real64), intent(out), target, dimension(:), optional :: work
3043 integer(int32), intent(out), optional :: olwork
3044 class(errors), intent(inout), optional, target :: err
3045
3046 ! Parameters
3047 real(real64), parameter :: zero = 0.0d0
3048
3049 ! Local Variables
3050 integer(int32) :: i, j, m, n, mn, k, istat, flag, lwork
3051 real(real64), pointer, dimension(:) :: wptr
3052 real(real64), allocatable, target, dimension(:) :: wrk
3053 real(real64), dimension(1) :: temp
3054 class(errors), pointer :: errmgr
3055 type(errors), target :: deferr
3056 character(len = :), allocatable :: errmsg
3057
3058 ! Initialization
3059 m = size(l, 1)
3060 n = size(l, 2)
3061 mn = min(m, n)
3062 if (present(err)) then
3063 errmgr => err
3064 else
3065 errmgr => deferr
3066 end if
3067
3068 ! Input Check
3069 flag = 0
3070 if (m > n) then
3071 flag = 1
3072 else if (size(tau) /= mn) then
3073 flag = 2
3074 else if (size(q, 1) /= n .or. size(q, 2) /= n) then
3075 flag = 3
3076 end if
3077 if (flag /= 0) then
3078 ! ERROR: One of the input arrays is not sized correctly
3079 allocate(character(len = 256) :: errmsg)
3080 write(errmsg, 100) "Input number ", flag, &
3081 " is not sized correctly."
3082 call errmgr%report_error("form_lq_no_pivot", trim(errmsg), &
3083 la_array_size_error)
3084 return
3085 end if
3086
3087 ! Workspace Query
3088 call dorglq(n, n, mn, q, n, tau, temp, -1, flag)
3089 lwork = int(temp(1), int32)
3090 if (present(olwork)) then
3091 olwork = lwork
3092 return
3093 end if
3094
3095 ! Local Memory Allocation
3096 if (present(work)) then
3097 if (size(work) < lwork) then
3098 ! ERROR: WORK not sized correctly
3099 call errmgr%report_error("form_lq_no_pivot", &
3100 "Incorrectly sized input array WORK, argument 4.", &
3101 la_array_size_error)
3102 return
3103 end if
3104 wptr => work(1:lwork)
3105 else
3106 allocate(wrk(lwork), stat = istat)
3107 if (istat /= 0) then
3108 ! ERROR: Out of memory
3109 call errmgr%report_error("form_lq_no_pivot", &
3110 "Insufficient memory available.", &
3111 la_out_of_memory_error)
3112 return
3113 end if
3114 wptr => wrk
3115 end if
3116
3117 ! Copy the upper triangular portion of L to Q, and then zero it out in L
3118 do j = 2, n
3119 k = min(j - 1, m)
3120 q(1:k,j) = l(1:k,j)
3121 l(1:k,j) = zero
3122 end do
3123
3124 ! Build Q
3125 call dorglq(n, n, mn, q, n, tau, wptr, lwork, flag)
3126
3127 ! Formatting
3128100 format(a, i0, a)
3129 end subroutine
3130
3131! ------------------------------------------------------------------------------
3132 module subroutine form_lq_no_pivot_cmplx(l, tau, q, work, olwork, err)
3133 ! Arguments
3134 complex(real64), intent(inout), dimension(:,:) :: l
3135 complex(real64), intent(in), dimension(:) :: tau
3136 complex(real64), intent(out), dimension(:,:) :: q
3137 complex(real64), intent(out), target, dimension(:), optional :: work
3138 integer(int32), intent(out), optional :: olwork
3139 class(errors), intent(inout), optional, target :: err
3140
3141 ! Parameters
3142 complex(real64), parameter :: zero = (0.0d0, 0.0d0)
3143
3144 ! Local Variables
3145 integer(int32) :: i, j, m, n, mn, k, istat, flag, lwork
3146 complex(real64), pointer, dimension(:) :: wptr
3147 complex(real64), allocatable, target, dimension(:) :: wrk
3148 complex(real64), dimension(1) :: temp
3149 class(errors), pointer :: errmgr
3150 type(errors), target :: deferr
3151 character(len = :), allocatable :: errmsg
3152
3153 ! Initialization
3154 m = size(l, 1)
3155 n = size(l, 2)
3156 mn = min(m, n)
3157 if (present(err)) then
3158 errmgr => err
3159 else
3160 errmgr => deferr
3161 end if
3162
3163 ! Input Check
3164 flag = 0
3165 if (m > n) then
3166 flag = 1
3167 else if (size(tau) /= mn) then
3168 flag = 2
3169 else if (size(q, 1) /= n .or. size(q, 2) /= n) then
3170 flag = 3
3171 end if
3172 if (flag /= 0) then
3173 ! ERROR: One of the input arrays is not sized correctly
3174 allocate(character(len = 256) :: errmsg)
3175 write(errmsg, 100) "Input number ", flag, &
3176 " is not sized correctly."
3177 call errmgr%report_error("form_lq_no_pivot_cmplx", trim(errmsg), &
3178 la_array_size_error)
3179 return
3180 end if
3181
3182 ! Workspace Query
3183 call zunglq(n, n, mn, q, n, tau, temp, -1, flag)
3184 lwork = int(temp(1), int32)
3185 if (present(olwork)) then
3186 olwork = lwork
3187 return
3188 end if
3189
3190 ! Local Memory Allocation
3191 if (present(work)) then
3192 if (size(work) < lwork) then
3193 ! ERROR: WORK not sized correctly
3194 call errmgr%report_error("form_lq_no_pivot_cmplx", &
3195 "Incorrectly sized input array WORK, argument 4.", &
3196 la_array_size_error)
3197 return
3198 end if
3199 wptr => work(1:lwork)
3200 else
3201 allocate(wrk(lwork), stat = istat)
3202 if (istat /= 0) then
3203 ! ERROR: Out of memory
3204 call errmgr%report_error("form_lq_no_pivot_cmplx", &
3205 "Insufficient memory available.", &
3206 la_out_of_memory_error)
3207 return
3208 end if
3209 wptr => wrk
3210 end if
3211
3212 ! Copy the upper triangular portion of L to Q, and then zero it out in L
3213 do j = 2, n
3214 k = min(j - 1, m)
3215 q(1:k,j) = l(1:k,j)
3216 l(1:k,j) = zero
3217 end do
3218
3219 ! Build Q
3220 call zunglq(n, n, mn, q, n, tau, wptr, lwork, flag)
3221
3222 ! Formatting
3223100 format(a, i0, a)
3224 end subroutine
3225
3226! ------------------------------------------------------------------------------
3227 module subroutine mult_lq_mtx(lside, trans, a, tau, c, work, olwork, err)
3228 ! Arguments
3229 logical, intent(in) :: lside, trans
3230 real(real64), intent(in), dimension(:,:) :: a
3231 real(real64), intent(in), dimension(:) :: tau
3232 real(real64), intent(inout), dimension(:,:) :: c
3233 real(real64), intent(out), target, dimension(:), optional :: work
3234 integer(int32), intent(out), optional :: olwork
3235 class(errors), intent(inout), optional, target :: err
3236
3237 ! Local Variables
3238 character :: side, t
3239 integer(int32) :: m, n, k, ncola, istat, flag, lwork
3240 real(real64), pointer, dimension(:) :: wptr
3241 real(real64), allocatable, target, dimension(:) :: wrk
3242 real(real64), dimension(1) :: temp
3243 class(errors), pointer :: errmgr
3244 type(errors), target :: deferr
3245 character(len = :), allocatable :: errmsg
3246
3247 ! Initialization
3248 m = size(c, 1)
3249 n = size(c, 2)
3250 k = size(tau)
3251 if (lside) then
3252 side = 'L'
3253 ncola = m
3254 else
3255 side = 'R'
3256 ncola = n
3257 end if
3258 if (trans) then
3259 t = 'T'
3260 else
3261 t = 'N'
3262 end if
3263 if (present(err)) then
3264 errmgr => err
3265 else
3266 errmgr => deferr
3267 end if
3268
3269 ! Input Check
3270 flag = 0
3271 if (size(a, 1) /= k .or. size(a, 2) /= ncola) then
3272 flag = 3
3273 end if
3274 if (flag /= 0) then
3275 ! ERROR: One of the input arrays is not sized correctly
3276 allocate(character(len = 256) :: errmsg)
3277 write(errmsg, 100) "Input number ", flag, &
3278 " is not sized correctly."
3279 call errmgr%report_error("mult_lq_mtx", trim(errmsg), &
3280 la_array_size_error)
3281 return
3282 end if
3283
3284 ! Workspace Query
3285 call dormlq(side, t, m, n, k, a, k, tau, c, m, temp, -1, flag)
3286 lwork = int(temp(1), int32)
3287 if (present(olwork)) then
3288 olwork = lwork
3289 return
3290 end if
3291
3292 ! Local Memory Allocation
3293 if (present(work)) then
3294 if (size(work) < lwork) then
3295 ! ERROR: WORK not sized correctly
3296 call errmgr%report_error("mult_lq_mtx", &
3297 "Incorrectly sized input array WORK, argument 6.", &
3298 la_array_size_error)
3299 return
3300 end if
3301 wptr => work(1:lwork)
3302 else
3303 allocate(wrk(lwork), stat = istat)
3304 if (istat /= 0) then
3305 ! ERROR: Out of memory
3306 call errmgr%report_error("mult_lq_mtx", &
3307 "Insufficient memory available.", &
3308 la_out_of_memory_error)
3309 return
3310 end if
3311 wptr => wrk
3312 end if
3313
3314 ! Call DORMLQ
3315 call dormlq(side, t, m, n, k, a, k, tau, c, m, wptr, lwork, flag)
3316
3317 ! Formatting
3318100 format(a, i0, a)
3319 end subroutine
3320
3321! ------------------------------------------------------------------------------
3322 module subroutine mult_lq_mtx_cmplx(lside, trans, a, tau, c, work, olwork, err)
3323 ! Arguments
3324 logical, intent(in) :: lside, trans
3325 complex(real64), intent(in), dimension(:,:) :: a
3326 complex(real64), intent(in), dimension(:) :: tau
3327 complex(real64), intent(inout), dimension(:,:) :: c
3328 complex(real64), intent(out), target, dimension(:), optional :: work
3329 integer(int32), intent(out), optional :: olwork
3330 class(errors), intent(inout), optional, target :: err
3331
3332 ! Local Variables
3333 character :: side, t
3334 integer(int32) :: m, n, k, ncola, istat, flag, lwork
3335 complex(real64), pointer, dimension(:) :: wptr
3336 complex(real64), allocatable, target, dimension(:) :: wrk
3337 complex(real64), dimension(1) :: temp
3338 class(errors), pointer :: errmgr
3339 type(errors), target :: deferr
3340 character(len = :), allocatable :: errmsg
3341
3342 ! Initialization
3343 m = size(c, 1)
3344 n = size(c, 2)
3345 k = size(tau)
3346 if (lside) then
3347 side = 'L'
3348 ncola = m
3349 else
3350 side = 'R'
3351 ncola = n
3352 end if
3353 if (trans) then
3354 t = 'C'
3355 else
3356 t = 'N'
3357 end if
3358 if (present(err)) then
3359 errmgr => err
3360 else
3361 errmgr => deferr
3362 end if
3363
3364 ! Input Check
3365 flag = 0
3366 if (size(a, 1) /= k .or. size(a, 2) /= ncola) then
3367 flag = 3
3368 end if
3369 if (flag /= 0) then
3370 ! ERROR: One of the input arrays is not sized correctly
3371 allocate(character(len = 256) :: errmsg)
3372 write(errmsg, 100) "Input number ", flag, &
3373 " is not sized correctly."
3374 call errmgr%report_error("mult_lq_mtx_cmplx", trim(errmsg), &
3375 la_array_size_error)
3376 return
3377 end if
3378
3379 ! Workspace Query
3380 call zunmlq(side, t, m, n, k, a, k, tau, c, m, temp, -1, flag)
3381 lwork = int(temp(1), int32)
3382 if (present(olwork)) then
3383 olwork = lwork
3384 return
3385 end if
3386
3387 ! Local Memory Allocation
3388 if (present(work)) then
3389 if (size(work) < lwork) then
3390 ! ERROR: WORK not sized correctly
3391 call errmgr%report_error("mult_lq_mtx_cmplx", &
3392 "Incorrectly sized input array WORK, argument 6.", &
3393 la_array_size_error)
3394 return
3395 end if
3396 wptr => work(1:lwork)
3397 else
3398 allocate(wrk(lwork), stat = istat)
3399 if (istat /= 0) then
3400 ! ERROR: Out of memory
3401 call errmgr%report_error("mult_lq_mtx_cmplx", &
3402 "Insufficient memory available.", &
3403 la_out_of_memory_error)
3404 return
3405 end if
3406 wptr => wrk
3407 end if
3408
3409 ! Call ZUNMLQ
3410 call zunmlq(side, t, m, n, k, a, k, tau, c, m, wptr, lwork, flag)
3411
3412 ! Formatting
3413100 format(a, i0, a)
3414 end subroutine
3415
3416! ------------------------------------------------------------------------------
3417 module subroutine mult_lq_vec(trans, a, tau, c, work, olwork, err)
3418 ! Arguments
3419 logical, intent(in) :: trans
3420 real(real64), intent(in), dimension(:,:) :: a
3421 real(real64), intent(in), dimension(:) :: tau
3422 real(real64), intent(inout), dimension(:) :: c
3423 real(real64), intent(out), target, dimension(:), optional :: work
3424 integer(int32), intent(out), optional :: olwork
3425 class(errors), intent(inout), optional, target :: err
3426
3427 ! Local Variables
3428 character :: side, t
3429 integer(int32) :: m, n, k, istat, flag, lwork
3430 real(real64), pointer, dimension(:) :: wptr
3431 real(real64), allocatable, target, dimension(:) :: wrk
3432 real(real64), dimension(1) :: temp
3433 class(errors), pointer :: errmgr
3434 type(errors), target :: deferr
3435 character(len = :), allocatable :: errmsg
3436
3437 ! Initialization
3438 m = size(c)
3439 n = 1
3440 k = size(tau)
3441 side = 'L'
3442 if (trans) then
3443 t = 'T'
3444 else
3445 t = 'N'
3446 end if
3447 if (present(err)) then
3448 errmgr => err
3449 else
3450 errmgr => deferr
3451 end if
3452
3453 ! Input Check
3454 flag = 0
3455 if (size(a, 1) /= k .or. size(a, 2) /= m) then
3456 flag = 3
3457 end if
3458 if (flag /= 0) then
3459 ! ERROR: One of the input arrays is not sized correctly
3460 allocate(character(len = 256) :: errmsg)
3461 write(errmsg, 100) "Input number ", flag, &
3462 " is not sized correctly."
3463 call errmgr%report_error("mult_lq_vec", trim(errmsg), &
3464 la_array_size_error)
3465 return
3466 end if
3467
3468 ! Workspace Query
3469 call dormlq(side, t, m, n, k, a, k, tau, c, m, temp, -1, flag)
3470 lwork = int(temp(1), int32)
3471 if (present(olwork)) then
3472 olwork = lwork
3473 return
3474 end if
3475
3476 ! Local Memory Allocation
3477 if (present(work)) then
3478 if (size(work) < lwork) then
3479 ! ERROR: WORK not sized correctly
3480 call errmgr%report_error("mult_lq_vec", &
3481 "Incorrectly sized input array WORK, argument 6.", &
3482 la_array_size_error)
3483 return
3484 end if
3485 wptr => work(1:lwork)
3486 else
3487 allocate(wrk(lwork), stat = istat)
3488 if (istat /= 0) then
3489 ! ERROR: Out of memory
3490 call errmgr%report_error("mult_lq_vec", &
3491 "Insufficient memory available.", &
3492 la_out_of_memory_error)
3493 return
3494 end if
3495 wptr => wrk
3496 end if
3497
3498 ! Call DORMLQ
3499 call dormlq(side, t, m, n, k, a, k, tau, c, m, wptr, lwork, flag)
3500
3501 ! Formatting
3502100 format(a, i0, a)
3503 end subroutine
3504
3505! ------------------------------------------------------------------------------
3506 module subroutine mult_lq_vec_cmplx(trans, a, tau, c, work, olwork, err)
3507 ! Arguments
3508 logical, intent(in) :: trans
3509 complex(real64), intent(in), dimension(:,:) :: a
3510 complex(real64), intent(in), dimension(:) :: tau
3511 complex(real64), intent(inout), dimension(:) :: c
3512 complex(real64), intent(out), target, dimension(:), optional :: work
3513 integer(int32), intent(out), optional :: olwork
3514 class(errors), intent(inout), optional, target :: err
3515
3516 ! Local Variables
3517 character :: side, t
3518 integer(int32) :: m, n, k, istat, flag, lwork
3519 complex(real64), pointer, dimension(:) :: wptr
3520 complex(real64), allocatable, target, dimension(:) :: wrk
3521 complex(real64), dimension(1) :: temp
3522 class(errors), pointer :: errmgr
3523 type(errors), target :: deferr
3524 character(len = :), allocatable :: errmsg
3525
3526 ! Initialization
3527 m = size(c)
3528 n = 1
3529 k = size(tau)
3530 side = 'L'
3531 if (trans) then
3532 t = 'C'
3533 else
3534 t = 'N'
3535 end if
3536 if (present(err)) then
3537 errmgr => err
3538 else
3539 errmgr => deferr
3540 end if
3541
3542 ! Input Check
3543 flag = 0
3544 if (size(a, 1) /= k .or. size(a, 2) /= m) then
3545 flag = 3
3546 end if
3547 if (flag /= 0) then
3548 ! ERROR: One of the input arrays is not sized correctly
3549 allocate(character(len = 256) :: errmsg)
3550 write(errmsg, 100) "Input number ", flag, &
3551 " is not sized correctly."
3552 call errmgr%report_error("mult_lq_vec_cmplx", trim(errmsg), &
3553 la_array_size_error)
3554 return
3555 end if
3556
3557 ! Workspace Query
3558 call zunmlq(side, t, m, n, k, a, k, tau, c, m, temp, -1, flag)
3559 lwork = int(temp(1), int32)
3560 if (present(olwork)) then
3561 olwork = lwork
3562 return
3563 end if
3564
3565 ! Local Memory Allocation
3566 if (present(work)) then
3567 if (size(work) < lwork) then
3568 ! ERROR: WORK not sized correctly
3569 call errmgr%report_error("mult_lq_vec_cmplx", &
3570 "Incorrectly sized input array WORK, argument 6.", &
3571 la_array_size_error)
3572 return
3573 end if
3574 wptr => work(1:lwork)
3575 else
3576 allocate(wrk(lwork), stat = istat)
3577 if (istat /= 0) then
3578 ! ERROR: Out of memory
3579 call errmgr%report_error("mult_lq_vec_cmplx", &
3580 "Insufficient memory available.", &
3581 la_out_of_memory_error)
3582 return
3583 end if
3584 wptr => wrk
3585 end if
3586
3587 ! Call ZUNMLQ
3588 call zunmlq(side, t, m, n, k, a, k, tau, c, m, wptr, lwork, flag)
3589
3590 ! Formatting
3591100 format(a, i0, a)
3592 end subroutine
3593
3594! ------------------------------------------------------------------------------
3595end submodule
linalg
Provides a set of common linear algebra routines.
Definition linalg.f90:145
qrupdate
A module providing explicit interfaces for the QRUPDATE library.
Definition qrupdate.f90:2