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_sorting.f90
1! linalg_sorting.f90
2
7submodule(linalg) linalg_sorting
8 use lapack
9 implicit none
10contains
11! ******************************************************************************
12! SORTING ROUTINES
13! ------------------------------------------------------------------------------
14 module subroutine sort_dbl_array(x, ascend)
15 ! Arguments
16 real(real64), intent(inout), dimension(:) :: x
17 logical, intent(in), optional :: ascend
18
19 ! Local Variables
20 character :: id
21 integer(int32) :: n, info
22
23 ! Initialization
24 if (present(ascend)) then
25 if (ascend) then
26 id = 'I'
27 else
28 id = 'D'
29 end if
30 else
31 id = 'I'
32 end if
33 n = size(x)
34
35 ! Process
36 call dlasrt(id, n, x, info)
37 end subroutine
38
39! ------------------------------------------------------------------------------
40 module subroutine sort_dbl_array_ind(x, ind, ascend, err)
41 ! Arguments
42 real(real64), intent(inout), dimension(:) :: x
43 integer(int32), intent(inout), dimension(:) :: ind
44 logical, intent(in), optional :: ascend
45 class(errors), intent(inout), optional, target :: err
46
47 ! Local Variables
48 class(errors), pointer :: errmgr
49 type(errors), target :: deferr
50 character(len = :), allocatable :: errmsg
51 integer(int32) :: n
52 logical :: dir
53
54 ! Initialization
55 n = size(x)
56 if (present(err)) then
57 errmgr => err
58 else
59 errmgr => deferr
60 end if
61 if (present(ascend)) then
62 dir = ascend
63 else
64 dir = .true. ! Ascend == true
65 end if
66
67 ! Input Check
68 if (size(ind) /= n) then
69 allocate(character(len = 256) :: errmsg)
70 write(errmsg, 100) &
71 "Expected the tracking array to be of size ", n, &
72 ", but found an array of size ", size(ind), "."
73 call errmgr%report_error("sort_dbl_array_ind", trim(errmsg), &
74 la_array_size_error)
75 return
76 end if
77 if (n <= 1) return
78
79 ! Process
80 call qsort_dbl_ind(dir, x, ind)
81
82 ! Formatting
83100 format(a, i0, a, i0, a)
84 end subroutine
85
86! ------------------------------------------------------------------------------
87 module subroutine sort_cmplx_array(x, ascend)
88 ! Arguments
89 complex(real64), intent(inout), dimension(:) :: x
90 logical, intent(in), optional :: ascend
91
92 ! Local Variables
93 logical :: dir
94
95 ! Initialization
96 if (present(ascend)) then
97 dir = ascend
98 else
99 dir = .true.
100 end if
101
102 ! Process
103 call qsort_cmplx(dir, x)
104 end subroutine
105
106! ------------------------------------------------------------------------------
107 module subroutine sort_cmplx_array_ind(x, ind, ascend, err)
108 ! Arguments
109 complex(real64), intent(inout), dimension(:) :: x
110 integer(int32), intent(inout), dimension(:) :: ind
111 logical, intent(in), optional :: ascend
112 class(errors), intent(inout), optional, target :: err
113
114 ! Local Variables
115 class(errors), pointer :: errmgr
116 type(errors), target :: deferr
117 character(len = :), allocatable :: errmsg
118 integer(int32) :: n
119 logical :: dir
120
121 ! Initialization
122 n = size(x)
123 if (present(err)) then
124 errmgr => err
125 else
126 errmgr => deferr
127 end if
128 if (present(ascend)) then
129 dir = ascend
130 else
131 dir = .true. ! Ascend == true
132 end if
133
134 ! Input Check
135 if (size(ind) /= n) then
136 allocate(character(len = 256) :: errmsg)
137 write(errmsg, 100) &
138 "Expected the tracking array to be of size ", n, &
139 ", but found an array of size ", size(ind), "."
140 call errmgr%report_error("sort_cmplx_array_ind", trim(errmsg), &
141 la_array_size_error)
142 return
143 end if
144 if (n <= 1) return
145
146 ! Process
147 call qsort_cmplx_ind(dir, x, ind)
148
149 ! Formatting
150100 format(a, i0, a, i0, a)
151 end subroutine
152
153! ------------------------------------------------------------------------------
154 module subroutine sort_eigen_cmplx(vals, vecs, ascend, err)
155 ! Arguments
156 complex(real64), intent(inout), dimension(:) :: vals
157 complex(real64), intent(inout), dimension(:,:) :: vecs
158 logical, intent(in), optional :: ascend
159 class(errors), intent(inout), optional, target :: err
160
161 ! Local Variables
162 class(errors), pointer :: errmgr
163 type(errors), target :: deferr
164 character(len = :), allocatable :: errmsg
165 integer(int32) :: i, n, flag
166 logical :: dir
167 integer(int32), allocatable, dimension(:) :: ind
168
169 ! Initialization
170 if (present(err)) then
171 errmgr => err
172 else
173 errmgr => deferr
174 end if
175 if (present(ascend)) then
176 dir = ascend
177 else
178 dir = .true. ! Ascend == true
179 end if
180
181 ! Ensure the eigenvector matrix is sized appropriately
182 n = size(vals)
183 if (size(vecs, 1) /= n .or. size(vecs, 2) /= n) then
184 ! ARRAY SIZE ERROR
185 allocate(character(len = 256) :: errmsg)
186 write(errmsg, 100) &
187 "Expected the eigenvector matrix to be of size ", n, &
188 "-by-", n, ", but found a matrix of size ", size(vecs, 1), &
189 "-by-", size(vecs, 2), "."
190 call errmgr%report_error("sort_eigen_cmplx", trim(errmsg), &
191 la_array_size_error)
192 end if
193
194 ! Allocate memory for the tracking array
195 allocate(ind(n), stat = flag)
196 if (flag /= 0) then
197 call errmgr%report_error("sort_eigen_cmplx", &
198 "Insufficient memory available.", la_out_of_memory_error)
199 return
200 end if
201 do i = 1, n
202 ind(i) = i
203 end do
204
205 ! Sort
206 call qsort_cmplx_ind(dir, vals, ind)
207
208 ! Shift the eigenvectors around to keep them associated with the
209 ! appropriate eigenvalue
210 vecs = vecs(:,ind)
211
212 ! Formatting
213100 format(a, i0, a, i0, a, i0, a, i0, a)
214 end subroutine
215
216! ------------------------------------------------------------------------------
217 module subroutine sort_eigen_dbl(vals, vecs, ascend, err)
218 ! Arguments
219 real(real64), intent(inout), dimension(:) :: vals
220 real(real64), intent(inout), dimension(:,:) :: vecs
221 logical, intent(in), optional :: ascend
222 class(errors), intent(inout), optional, target :: err
223
224 ! Local Variables
225 class(errors), pointer :: errmgr
226 type(errors), target :: deferr
227 character(len = :), allocatable :: errmsg
228 integer(int32) :: i, n, flag
229 logical :: dir
230 integer(int32), allocatable, dimension(:) :: ind
231
232 ! Initialization
233 if (present(err)) then
234 errmgr => err
235 else
236 errmgr => deferr
237 end if
238 if (present(ascend)) then
239 dir = ascend
240 else
241 dir = .true. ! Ascend == true
242 end if
243
244 ! Ensure the eigenvector matrix is sized appropriately
245 n = size(vals)
246 if (size(vecs, 1) /= n .or. size(vecs, 2) /= n) then
247 ! ARRAY SIZE ERROR
248 allocate(character(len = 256) :: errmsg)
249 write(errmsg, 100) &
250 "Expected the eigenvector matrix to be of size ", n, &
251 "-by-", n, ", but found a matrix of size ", size(vecs, 1), &
252 "-by-", size(vecs, 2), "."
253 call errmgr%report_error("sort_eigen_dbl", trim(errmsg), &
254 la_array_size_error)
255 end if
256
257 ! Allocate memory for the tracking array
258 allocate(ind(n), stat = flag)
259 if (flag /= 0) then
260 call errmgr%report_error("sort_eigen_dbl", &
261 "Insufficient memory available.", la_out_of_memory_error)
262 return
263 end if
264 do i = 1, n
265 ind(i) = i
266 end do
267
268 ! Sort
269 call qsort_dbl_ind(dir, vals, ind)
270
271 ! Shift the eigenvectors around to keep them associated with the
272 ! appropriate eigenvalue
273 vecs = vecs(:,ind)
274
275 ! Formatting
276100 format(a, i0, a, i0, a, i0, a, i0, a)
277 end subroutine
278
279! ------------------------------------------------------------------------------
280 module subroutine sort_int32_array(x, ascend)
281 ! Arguments
282 integer(int32), intent(inout), dimension(:) :: x
283 logical, intent(in), optional :: ascend
284
285 ! Local Variables
286 logical :: dir
287
288 ! Initialization
289 if (present(ascend)) then
290 dir = ascend
291 else
292 dir = .true.
293 end if
294
295 ! Process
296 call qsort_int32(dir, x)
297 end subroutine
298
299! ------------------------------------------------------------------------------
300 module subroutine sort_int32_array_ind(x, ind, ascend, err)
301 ! Arguments
302 integer(int32), intent(inout), dimension(:) :: x
303 integer(int32), intent(inout), dimension(:) :: ind
304 logical, intent(in), optional :: ascend
305 class(errors), intent(inout), optional, target :: err
306
307 ! Local Variables
308 class(errors), pointer :: errmgr
309 type(errors), target :: deferr
310 character(len = :), allocatable :: errmsg
311 integer(int32) :: n
312 logical :: dir
313
314 ! Initialization
315 n = size(x)
316 if (present(err)) then
317 errmgr => err
318 else
319 errmgr => deferr
320 end if
321 if (present(ascend)) then
322 dir = ascend
323 else
324 dir = .true. ! Ascend == true
325 end if
326
327 ! Input Check
328 if (size(ind) /= n) then
329 allocate(character(len = 256) :: errmsg)
330 write(errmsg, 100) &
331 "Expected the tracking array to be of size ", n, &
332 ", but found an array of size ", size(ind), "."
333 call errmgr%report_error("sort_int32_array_ind", trim(errmsg), &
334 la_array_size_error)
335 return
336 end if
337 if (n <= 1) return
338
339 ! Process
340 call qsort_int32_ind(dir, x, ind)
341
342 ! Formatting
343100 format(a, i0, a, i0, a)
344 end subroutine
345
346! ******************************************************************************
347! PRIVATE HELPER ROUTINES
348! ------------------------------------------------------------------------------
362 recursive subroutine qsort_dbl_ind(ascend, x, ind)
363 ! Arguments
364 logical, intent(in) :: ascend
365 real(real64), intent(inout), dimension(:) :: x
366 integer(int32), intent(inout), dimension(:) :: ind
367
368 ! Local Variables
369 integer(int32) :: iq
370
371 ! Process
372 if (size(x) > 1) then
373 call dbl_partition_ind(ascend, x, ind, iq)
374 call qsort_dbl_ind(ascend, x(:iq-1), ind(:iq-1))
375 call qsort_dbl_ind(ascend, x(iq:), ind(iq:))
376 end if
377 end subroutine
378
379! ------------------------------------------------------------------------------
395 subroutine dbl_partition_ind(ascend, x, ind, marker)
396 ! Arguments
397 logical, intent(in) :: ascend
398 real(real64), intent(inout), dimension(:) :: x
399 integer(int32), intent(inout), dimension(:) :: ind
400 integer(int32), intent(out) :: marker
401
402 ! Local Variables
403 integer(int32) :: i, j, itemp
404 real(real64) :: temp, pivot
405
406 ! Process
407 pivot = x(1)
408 i = 0
409 j = size(x) + 1
410 if (ascend) then
411 ! Ascending Sort
412 do
413 j = j - 1
414 do
415 if (x(j) <= pivot) exit
416 j = j - 1
417 end do
418 i = i + 1
419 do
420 if (x(i) >= pivot) exit
421 i = i + 1
422 end do
423 if (i < j) then
424 ! Exchage X(I) and X(J)
425 temp = x(i)
426 x(i) = x(j)
427 x(j) = temp
428
429 itemp = ind(i)
430 ind(i) = ind(j)
431 ind(j) = itemp
432 else if (i == j) then
433 marker = i + 1
434 return
435 else
436 marker = i
437 return
438 end if
439 end do
440 else
441 ! Descending Sort
442 do
443 j = j - 1
444 do
445 if (x(j) >= pivot) exit
446 j = j - 1
447 end do
448 i = i + 1
449 do
450 if (x(i) <= pivot) exit
451 i = i + 1
452 end do
453 if (i < j) then
454 ! Exchage X(I) and X(J)
455 temp = x(i)
456 x(i) = x(j)
457 x(j) = temp
458
459 itemp = ind(i)
460 ind(i) = ind(j)
461 ind(j) = itemp
462 else if (i == j) then
463 marker = i + 1
464 return
465 else
466 marker = i
467 return
468 end if
469 end do
470 end if
471 end subroutine
472
473! ------------------------------------------------------------------------------
488 recursive subroutine qsort_cmplx(ascend, x)
489 ! Arguments
490 logical, intent(in) :: ascend
491 complex(real64), intent(inout), dimension(:) :: x
492
493 ! Local Variables
494 integer(int32) :: iq
495
496 ! Process
497 if (size(x) > 1) then
498 call cmplx_partition(ascend, x, iq)
499 call qsort_cmplx(ascend, x(:iq-1))
500 call qsort_cmplx(ascend, x(iq:))
501 end if
502 end subroutine
503
504! ------------------------------------------------------------------------------
521 subroutine cmplx_partition(ascend, x, marker)
522 ! Arguments
523 logical, intent(in) :: ascend
524 complex(real64), intent(inout), dimension(:) :: x
525 integer(int32), intent(out) :: marker
526
527 ! Local Variables
528 integer(int32) :: i, j
529 complex(real64) :: temp
530 real(real64) :: pivot
531
532 ! Process
533 pivot = real(x(1), real64)
534 i = 0
535 j = size(x) + 1
536 if (ascend) then
537 ! Ascending Sort
538 do
539 j = j - 1
540 do
541 if (real(x(j), real64) <= pivot) exit
542 j = j - 1
543 end do
544 i = i + 1
545 do
546 if (real(x(i), real64) >= pivot) exit
547 i = i + 1
548 end do
549 if (i < j) then
550 ! Exchage X(I) and X(J)
551 temp = x(i)
552 x(i) = x(j)
553 x(j) = temp
554 else if (i == j) then
555 marker = i + 1
556 return
557 else
558 marker = i
559 return
560 end if
561 end do
562 else
563 ! Descending Sort
564 do
565 j = j - 1
566 do
567 if (real(x(j), real64) >= pivot) exit
568 j = j - 1
569 end do
570 i = i + 1
571 do
572 if (real(x(i), real64) <= pivot) exit
573 i = i + 1
574 end do
575 if (i < j) then
576 ! Exchage X(I) and X(J)
577 temp = x(i)
578 x(i) = x(j)
579 x(j) = temp
580 else if (i == j) then
581 marker = i + 1
582 return
583 else
584 marker = i
585 return
586 end if
587 end do
588 end if
589 end subroutine
590
591! ------------------------------------------------------------------------------
609 recursive subroutine qsort_cmplx_ind(ascend, x, ind)
610 ! Arguments
611 logical, intent(in) :: ascend
612 complex(real64), intent(inout), dimension(:) :: x
613 integer(int32), intent(inout), dimension(:) :: ind
614
615 ! Local Variables
616 integer(int32) :: iq
617
618 ! Process
619 if (size(x) > 1) then
620 call cmplx_partition_ind(ascend, x, ind, iq)
621 call qsort_cmplx_ind(ascend, x(:iq-1), ind(:iq-1))
622 call qsort_cmplx_ind(ascend, x(iq:), ind(iq:))
623 end if
624 end subroutine
625
626! ------------------------------------------------------------------------------
646 subroutine cmplx_partition_ind(ascend, x, ind, marker)
647 ! Arguments
648 logical, intent(in) :: ascend
649 complex(real64), intent(inout), dimension(:) :: x
650 integer(int32), intent(inout), dimension(:) :: ind
651 integer(int32), intent(out) :: marker
652
653 ! Local Variables
654 integer(int32) :: i, j, itemp
655 complex(real64) :: temp
656 real(real64) :: pivot
657
658 ! Process
659 pivot = real(x(1), real64)
660 i = 0
661 j = size(x) + 1
662 if (ascend) then
663 ! Ascending Sort
664 do
665 j = j - 1
666 do
667 if (real(x(j), real64) <= pivot) exit
668 j = j - 1
669 end do
670 i = i + 1
671 do
672 if (real(x(i), real64) >= pivot) exit
673 i = i + 1
674 end do
675 if (i < j) then
676 ! Exchage X(I) and X(J)
677 temp = x(i)
678 x(i) = x(j)
679 x(j) = temp
680
681 itemp = ind(i)
682 ind(i) = ind(j)
683 ind(j) = itemp
684 else if (i == j) then
685 marker = i + 1
686 return
687 else
688 marker = i
689 return
690 end if
691 end do
692 else
693 ! Descending Sort
694 do
695 j = j - 1
696 do
697 if (real(x(j), real64) >= pivot) exit
698 j = j - 1
699 end do
700 i = i + 1
701 do
702 if (real(x(i), real64) <= pivot) exit
703 i = i + 1
704 end do
705 if (i < j) then
706 ! Exchage X(I) and X(J)
707 temp = x(i)
708 x(i) = x(j)
709 x(j) = temp
710
711 itemp = ind(i)
712 ind(i) = ind(j)
713 ind(j) = itemp
714 else if (i == j) then
715 marker = i + 1
716 return
717 else
718 marker = i
719 return
720 end if
721 end do
722 end if
723 end subroutine
724
725! ------------------------------------------------------------------------------
736 recursive subroutine qsort_int32(ascend, x)
737 ! Arguments
738 logical, intent(in) :: ascend
739 integer(int32), intent(inout), dimension(:) :: x
740
741 ! Local Variables
742 integer(int32) :: iq
743
744 ! Process
745 if (size(x) > 1) then
746 call int32_partition(ascend, x, iq)
747 call qsort_int32(ascend, x(:iq-1))
748 call qsort_int32(ascend, x(iq:))
749 end if
750 end subroutine
751
752! ------------------------------------------------------------------------------
765 subroutine int32_partition(ascend, x, marker)
766 ! Arguments
767 logical, intent(in) :: ascend
768 integer(int32), intent(inout), dimension(:) :: x
769 integer(int32), intent(out) :: marker
770
771 ! Local Variables
772 integer(int32) :: i, j, temp, pivot
773
774 ! Process
775 pivot = x(1)
776 i = 0
777 j = size(x) + 1
778 if (ascend) then
779 ! Ascending Sort
780 do
781 j = j - 1
782 do
783 if (x(j) <= pivot) exit
784 j = j - 1
785 end do
786 i = i + 1
787 do
788 if (x(i) >= pivot) exit
789 i = i + 1
790 end do
791 if (i < j) then
792 ! Exchage X(I) and X(J)
793 temp = x(i)
794 x(i) = x(j)
795 x(j) = temp
796 else if (i == j) then
797 marker = i + 1
798 return
799 else
800 marker = i
801 return
802 end if
803 end do
804 else
805 ! Descending Sort
806 do
807 j = j - 1
808 do
809 if (x(j) >= pivot) exit
810 j = j - 1
811 end do
812 i = i + 1
813 do
814 if (x(i) <= pivot) exit
815 i = i + 1
816 end do
817 if (i < j) then
818 ! Exchage X(I) and X(J)
819 temp = x(i)
820 x(i) = x(j)
821 x(j) = temp
822 else if (i == j) then
823 marker = i + 1
824 return
825 else
826 marker = i
827 return
828 end if
829 end do
830 end if
831 end subroutine
832
833! ------------------------------------------------------------------------------
847 recursive subroutine qsort_int32_ind(ascend, x, ind)
848 ! Arguments
849 logical, intent(in) :: ascend
850 integer(int32), intent(inout), dimension(:) :: x
851 integer(int32), intent(inout), dimension(:) :: ind
852
853 ! Local Variables
854 integer(int32) :: iq
855
856 ! Process
857 if (size(x) > 1) then
858 call int32_partition_ind(ascend, x, ind, iq)
859 call qsort_int32_ind(ascend, x(:iq-1), ind(:iq-1))
860 call qsort_int32_ind(ascend, x(iq:), ind(iq:))
861 end if
862 end subroutine
863
864! ------------------------------------------------------------------------------
880 subroutine int32_partition_ind(ascend, x, ind, marker)
881 ! Arguments
882 logical, intent(in) :: ascend
883 integer(int32), intent(inout), dimension(:) :: x
884 integer(int32), intent(inout), dimension(:) :: ind
885 integer(int32), intent(out) :: marker
886
887 ! Local Variables
888 integer(int32) :: i, j, itemp, temp, pivot
889
890 ! Process
891 pivot = x(1)
892 i = 0
893 j = size(x) + 1
894 if (ascend) then
895 ! Ascending Sort
896 do
897 j = j - 1
898 do
899 if (x(j) <= pivot) exit
900 j = j - 1
901 end do
902 i = i + 1
903 do
904 if (x(i) >= pivot) exit
905 i = i + 1
906 end do
907 if (i < j) then
908 ! Exchage X(I) and X(J)
909 temp = x(i)
910 x(i) = x(j)
911 x(j) = temp
912
913 itemp = ind(i)
914 ind(i) = ind(j)
915 ind(j) = itemp
916 else if (i == j) then
917 marker = i + 1
918 return
919 else
920 marker = i
921 return
922 end if
923 end do
924 else
925 ! Descending Sort
926 do
927 j = j - 1
928 do
929 if (x(j) >= pivot) exit
930 j = j - 1
931 end do
932 i = i + 1
933 do
934 if (x(i) <= pivot) exit
935 i = i + 1
936 end do
937 if (i < j) then
938 ! Exchage X(I) and X(J)
939 temp = x(i)
940 x(i) = x(j)
941 x(j) = temp
942
943 itemp = ind(i)
944 ind(i) = ind(j)
945 ind(j) = itemp
946 else if (i == j) then
947 marker = i + 1
948 return
949 else
950 marker = i
951 return
952 end if
953 end do
954 end if
955 end subroutine
956
957! ------------------------------------------------------------------------------
958end submodule
linalg
Provides a set of common linear algebra routines.
Definition linalg.f90:145