mirror of
https://github.com/QuantumPackage/qp2.git
synced 2024-11-09 06:53:38 +01:00
This commit is contained in:
parent
d5fb21fe12
commit
050ee3af5d
@ -69,9 +69,9 @@ END_DOC
|
||||
if ( (scf_algorithm == 'DIIS').and.(dabs(Delta_energy_SCF) > 1.d-6) ) then
|
||||
|
||||
! Store Fock and error matrices at each iteration
|
||||
index_dim_DIIS = mod(dim_DIIS-1,max_dim_DIIS)+1
|
||||
do j=1,ao_num
|
||||
do i=1,ao_num
|
||||
index_dim_DIIS = mod(dim_DIIS-1,max_dim_DIIS)+1
|
||||
Fock_matrix_DIIS (i,j,index_dim_DIIS) = Fock_matrix_AO(i,j)
|
||||
error_matrix_DIIS(i,j,index_dim_DIIS) = FPS_SPF_matrix_AO(i,j)
|
||||
enddo
|
||||
|
||||
373
src/utils/qsort.c
Normal file
373
src/utils/qsort.c
Normal file
@ -0,0 +1,373 @@
|
||||
/* [[file:~/qp2/src/utils/qsort.org::*Generated%20C%20file][Generated C file:1]] */
|
||||
#include <stdlib.h>
|
||||
#include <stdint.h>
|
||||
|
||||
struct int16_t_comp {
|
||||
int16_t x;
|
||||
int32_t i;
|
||||
};
|
||||
|
||||
int compare_int16_t( const void * l, const void * r )
|
||||
{
|
||||
const int16_t * restrict _l= l;
|
||||
const int16_t * restrict _r= r;
|
||||
if( *_l > *_r ) return 1;
|
||||
if( *_l < *_r ) return -1;
|
||||
return 0;
|
||||
}
|
||||
|
||||
void qsort_int16_t(int16_t* restrict A_in, int32_t* restrict iorder, int32_t isize) {
|
||||
struct int16_t_comp* A = malloc(isize * sizeof(struct int16_t_comp));
|
||||
if (A == NULL) return;
|
||||
|
||||
for (int i=0 ; i<isize ; ++i) {
|
||||
A[i].x = A_in[i];
|
||||
A[i].i = iorder[i];
|
||||
}
|
||||
|
||||
qsort( (void*) A, (size_t) isize, sizeof(struct int16_t_comp), compare_int16_t);
|
||||
|
||||
for (int i=0 ; i<isize ; ++i) {
|
||||
A_in[i] = A[i].x;
|
||||
iorder[i] = A[i].i;
|
||||
}
|
||||
free(A);
|
||||
}
|
||||
|
||||
void qsort_int16_t_noidx(int16_t* A, int32_t isize) {
|
||||
qsort( (void*) A, (size_t) isize, sizeof(int16_t), compare_int16_t);
|
||||
}
|
||||
|
||||
|
||||
struct int16_t_comp_big {
|
||||
int16_t x;
|
||||
int64_t i;
|
||||
};
|
||||
|
||||
int compare_int16_t_big( const void * l, const void * r )
|
||||
{
|
||||
const int16_t * restrict _l= l;
|
||||
const int16_t * restrict _r= r;
|
||||
if( *_l > *_r ) return 1;
|
||||
if( *_l < *_r ) return -1;
|
||||
return 0;
|
||||
}
|
||||
|
||||
void qsort_int16_t_big(int16_t* restrict A_in, int64_t* restrict iorder, int64_t isize) {
|
||||
struct int16_t_comp_big* A = malloc(isize * sizeof(struct int16_t_comp_big));
|
||||
if (A == NULL) return;
|
||||
|
||||
for (int i=0 ; i<isize ; ++i) {
|
||||
A[i].x = A_in[i];
|
||||
A[i].i = iorder[i];
|
||||
}
|
||||
|
||||
qsort( (void*) A, (size_t) isize, sizeof(struct int16_t_comp_big), compare_int16_t_big);
|
||||
|
||||
for (int i=0 ; i<isize ; ++i) {
|
||||
A_in[i] = A[i].x;
|
||||
iorder[i] = A[i].i;
|
||||
}
|
||||
free(A);
|
||||
}
|
||||
|
||||
void qsort_int16_t_noidx_big(int16_t* A, int64_t isize) {
|
||||
qsort( (void*) A, (size_t) isize, sizeof(int16_t), compare_int16_t_big);
|
||||
}
|
||||
|
||||
|
||||
struct int32_t_comp {
|
||||
int32_t x;
|
||||
int32_t i;
|
||||
};
|
||||
|
||||
int compare_int32_t( const void * l, const void * r )
|
||||
{
|
||||
const int32_t * restrict _l= l;
|
||||
const int32_t * restrict _r= r;
|
||||
if( *_l > *_r ) return 1;
|
||||
if( *_l < *_r ) return -1;
|
||||
return 0;
|
||||
}
|
||||
|
||||
void qsort_int32_t(int32_t* restrict A_in, int32_t* restrict iorder, int32_t isize) {
|
||||
struct int32_t_comp* A = malloc(isize * sizeof(struct int32_t_comp));
|
||||
if (A == NULL) return;
|
||||
|
||||
for (int i=0 ; i<isize ; ++i) {
|
||||
A[i].x = A_in[i];
|
||||
A[i].i = iorder[i];
|
||||
}
|
||||
|
||||
qsort( (void*) A, (size_t) isize, sizeof(struct int32_t_comp), compare_int32_t);
|
||||
|
||||
for (int i=0 ; i<isize ; ++i) {
|
||||
A_in[i] = A[i].x;
|
||||
iorder[i] = A[i].i;
|
||||
}
|
||||
free(A);
|
||||
}
|
||||
|
||||
void qsort_int32_t_noidx(int32_t* A, int32_t isize) {
|
||||
qsort( (void*) A, (size_t) isize, sizeof(int32_t), compare_int32_t);
|
||||
}
|
||||
|
||||
|
||||
struct int32_t_comp_big {
|
||||
int32_t x;
|
||||
int64_t i;
|
||||
};
|
||||
|
||||
int compare_int32_t_big( const void * l, const void * r )
|
||||
{
|
||||
const int32_t * restrict _l= l;
|
||||
const int32_t * restrict _r= r;
|
||||
if( *_l > *_r ) return 1;
|
||||
if( *_l < *_r ) return -1;
|
||||
return 0;
|
||||
}
|
||||
|
||||
void qsort_int32_t_big(int32_t* restrict A_in, int64_t* restrict iorder, int64_t isize) {
|
||||
struct int32_t_comp_big* A = malloc(isize * sizeof(struct int32_t_comp_big));
|
||||
if (A == NULL) return;
|
||||
|
||||
for (int i=0 ; i<isize ; ++i) {
|
||||
A[i].x = A_in[i];
|
||||
A[i].i = iorder[i];
|
||||
}
|
||||
|
||||
qsort( (void*) A, (size_t) isize, sizeof(struct int32_t_comp_big), compare_int32_t_big);
|
||||
|
||||
for (int i=0 ; i<isize ; ++i) {
|
||||
A_in[i] = A[i].x;
|
||||
iorder[i] = A[i].i;
|
||||
}
|
||||
free(A);
|
||||
}
|
||||
|
||||
void qsort_int32_t_noidx_big(int32_t* A, int64_t isize) {
|
||||
qsort( (void*) A, (size_t) isize, sizeof(int32_t), compare_int32_t_big);
|
||||
}
|
||||
|
||||
|
||||
struct int64_t_comp {
|
||||
int64_t x;
|
||||
int32_t i;
|
||||
};
|
||||
|
||||
int compare_int64_t( const void * l, const void * r )
|
||||
{
|
||||
const int64_t * restrict _l= l;
|
||||
const int64_t * restrict _r= r;
|
||||
if( *_l > *_r ) return 1;
|
||||
if( *_l < *_r ) return -1;
|
||||
return 0;
|
||||
}
|
||||
|
||||
void qsort_int64_t(int64_t* restrict A_in, int32_t* restrict iorder, int32_t isize) {
|
||||
struct int64_t_comp* A = malloc(isize * sizeof(struct int64_t_comp));
|
||||
if (A == NULL) return;
|
||||
|
||||
for (int i=0 ; i<isize ; ++i) {
|
||||
A[i].x = A_in[i];
|
||||
A[i].i = iorder[i];
|
||||
}
|
||||
|
||||
qsort( (void*) A, (size_t) isize, sizeof(struct int64_t_comp), compare_int64_t);
|
||||
|
||||
for (int i=0 ; i<isize ; ++i) {
|
||||
A_in[i] = A[i].x;
|
||||
iorder[i] = A[i].i;
|
||||
}
|
||||
free(A);
|
||||
}
|
||||
|
||||
void qsort_int64_t_noidx(int64_t* A, int32_t isize) {
|
||||
qsort( (void*) A, (size_t) isize, sizeof(int64_t), compare_int64_t);
|
||||
}
|
||||
|
||||
|
||||
struct int64_t_comp_big {
|
||||
int64_t x;
|
||||
int64_t i;
|
||||
};
|
||||
|
||||
int compare_int64_t_big( const void * l, const void * r )
|
||||
{
|
||||
const int64_t * restrict _l= l;
|
||||
const int64_t * restrict _r= r;
|
||||
if( *_l > *_r ) return 1;
|
||||
if( *_l < *_r ) return -1;
|
||||
return 0;
|
||||
}
|
||||
|
||||
void qsort_int64_t_big(int64_t* restrict A_in, int64_t* restrict iorder, int64_t isize) {
|
||||
struct int64_t_comp_big* A = malloc(isize * sizeof(struct int64_t_comp_big));
|
||||
if (A == NULL) return;
|
||||
|
||||
for (int i=0 ; i<isize ; ++i) {
|
||||
A[i].x = A_in[i];
|
||||
A[i].i = iorder[i];
|
||||
}
|
||||
|
||||
qsort( (void*) A, (size_t) isize, sizeof(struct int64_t_comp_big), compare_int64_t_big);
|
||||
|
||||
for (int i=0 ; i<isize ; ++i) {
|
||||
A_in[i] = A[i].x;
|
||||
iorder[i] = A[i].i;
|
||||
}
|
||||
free(A);
|
||||
}
|
||||
|
||||
void qsort_int64_t_noidx_big(int64_t* A, int64_t isize) {
|
||||
qsort( (void*) A, (size_t) isize, sizeof(int64_t), compare_int64_t_big);
|
||||
}
|
||||
|
||||
|
||||
struct double_comp {
|
||||
double x;
|
||||
int32_t i;
|
||||
};
|
||||
|
||||
int compare_double( const void * l, const void * r )
|
||||
{
|
||||
const double * restrict _l= l;
|
||||
const double * restrict _r= r;
|
||||
if( *_l > *_r ) return 1;
|
||||
if( *_l < *_r ) return -1;
|
||||
return 0;
|
||||
}
|
||||
|
||||
void qsort_double(double* restrict A_in, int32_t* restrict iorder, int32_t isize) {
|
||||
struct double_comp* A = malloc(isize * sizeof(struct double_comp));
|
||||
if (A == NULL) return;
|
||||
|
||||
for (int i=0 ; i<isize ; ++i) {
|
||||
A[i].x = A_in[i];
|
||||
A[i].i = iorder[i];
|
||||
}
|
||||
|
||||
qsort( (void*) A, (size_t) isize, sizeof(struct double_comp), compare_double);
|
||||
|
||||
for (int i=0 ; i<isize ; ++i) {
|
||||
A_in[i] = A[i].x;
|
||||
iorder[i] = A[i].i;
|
||||
}
|
||||
free(A);
|
||||
}
|
||||
|
||||
void qsort_double_noidx(double* A, int32_t isize) {
|
||||
qsort( (void*) A, (size_t) isize, sizeof(double), compare_double);
|
||||
}
|
||||
|
||||
|
||||
struct double_comp_big {
|
||||
double x;
|
||||
int64_t i;
|
||||
};
|
||||
|
||||
int compare_double_big( const void * l, const void * r )
|
||||
{
|
||||
const double * restrict _l= l;
|
||||
const double * restrict _r= r;
|
||||
if( *_l > *_r ) return 1;
|
||||
if( *_l < *_r ) return -1;
|
||||
return 0;
|
||||
}
|
||||
|
||||
void qsort_double_big(double* restrict A_in, int64_t* restrict iorder, int64_t isize) {
|
||||
struct double_comp_big* A = malloc(isize * sizeof(struct double_comp_big));
|
||||
if (A == NULL) return;
|
||||
|
||||
for (int i=0 ; i<isize ; ++i) {
|
||||
A[i].x = A_in[i];
|
||||
A[i].i = iorder[i];
|
||||
}
|
||||
|
||||
qsort( (void*) A, (size_t) isize, sizeof(struct double_comp_big), compare_double_big);
|
||||
|
||||
for (int i=0 ; i<isize ; ++i) {
|
||||
A_in[i] = A[i].x;
|
||||
iorder[i] = A[i].i;
|
||||
}
|
||||
free(A);
|
||||
}
|
||||
|
||||
void qsort_double_noidx_big(double* A, int64_t isize) {
|
||||
qsort( (void*) A, (size_t) isize, sizeof(double), compare_double_big);
|
||||
}
|
||||
|
||||
|
||||
struct float_comp {
|
||||
float x;
|
||||
int32_t i;
|
||||
};
|
||||
|
||||
int compare_float( const void * l, const void * r )
|
||||
{
|
||||
const float * restrict _l= l;
|
||||
const float * restrict _r= r;
|
||||
if( *_l > *_r ) return 1;
|
||||
if( *_l < *_r ) return -1;
|
||||
return 0;
|
||||
}
|
||||
|
||||
void qsort_float(float* restrict A_in, int32_t* restrict iorder, int32_t isize) {
|
||||
struct float_comp* A = malloc(isize * sizeof(struct float_comp));
|
||||
if (A == NULL) return;
|
||||
|
||||
for (int i=0 ; i<isize ; ++i) {
|
||||
A[i].x = A_in[i];
|
||||
A[i].i = iorder[i];
|
||||
}
|
||||
|
||||
qsort( (void*) A, (size_t) isize, sizeof(struct float_comp), compare_float);
|
||||
|
||||
for (int i=0 ; i<isize ; ++i) {
|
||||
A_in[i] = A[i].x;
|
||||
iorder[i] = A[i].i;
|
||||
}
|
||||
free(A);
|
||||
}
|
||||
|
||||
void qsort_float_noidx(float* A, int32_t isize) {
|
||||
qsort( (void*) A, (size_t) isize, sizeof(float), compare_float);
|
||||
}
|
||||
|
||||
|
||||
struct float_comp_big {
|
||||
float x;
|
||||
int64_t i;
|
||||
};
|
||||
|
||||
int compare_float_big( const void * l, const void * r )
|
||||
{
|
||||
const float * restrict _l= l;
|
||||
const float * restrict _r= r;
|
||||
if( *_l > *_r ) return 1;
|
||||
if( *_l < *_r ) return -1;
|
||||
return 0;
|
||||
}
|
||||
|
||||
void qsort_float_big(float* restrict A_in, int64_t* restrict iorder, int64_t isize) {
|
||||
struct float_comp_big* A = malloc(isize * sizeof(struct float_comp_big));
|
||||
if (A == NULL) return;
|
||||
|
||||
for (int i=0 ; i<isize ; ++i) {
|
||||
A[i].x = A_in[i];
|
||||
A[i].i = iorder[i];
|
||||
}
|
||||
|
||||
qsort( (void*) A, (size_t) isize, sizeof(struct float_comp_big), compare_float_big);
|
||||
|
||||
for (int i=0 ; i<isize ; ++i) {
|
||||
A_in[i] = A[i].x;
|
||||
iorder[i] = A[i].i;
|
||||
}
|
||||
free(A);
|
||||
}
|
||||
|
||||
void qsort_float_noidx_big(float* A, int64_t isize) {
|
||||
qsort( (void*) A, (size_t) isize, sizeof(float), compare_float_big);
|
||||
}
|
||||
/* Generated C file:1 ends here */
|
||||
169
src/utils/qsort.org
Normal file
169
src/utils/qsort.org
Normal file
@ -0,0 +1,169 @@
|
||||
#+TITLE: Quick sort binding for Fortran
|
||||
|
||||
* C template
|
||||
|
||||
#+NAME: c_template
|
||||
#+BEGIN_SRC c
|
||||
struct TYPE_comp_big {
|
||||
TYPE x;
|
||||
int32_t i;
|
||||
};
|
||||
|
||||
int compare_TYPE_big( const void * l, const void * r )
|
||||
{
|
||||
const TYPE * restrict _l= l;
|
||||
const TYPE * restrict _r= r;
|
||||
if( *_l > *_r ) return 1;
|
||||
if( *_l < *_r ) return -1;
|
||||
return 0;
|
||||
}
|
||||
|
||||
void qsort_TYPE_big(TYPE* restrict A_in, int32_t* restrict iorder, int32_t isize) {
|
||||
struct TYPE_comp_big* A = malloc(isize * sizeof(struct TYPE_comp_big));
|
||||
if (A == NULL) return;
|
||||
|
||||
for (int i=0 ; i<isize ; ++i) {
|
||||
A[i].x = A_in[i];
|
||||
A[i].i = iorder[i];
|
||||
}
|
||||
|
||||
qsort( (void*) A, (size_t) isize, sizeof(struct TYPE_comp_big), compare_TYPE_big);
|
||||
|
||||
for (int i=0 ; i<isize ; ++i) {
|
||||
A_in[i] = A[i].x;
|
||||
iorder[i] = A[i].i;
|
||||
}
|
||||
free(A);
|
||||
}
|
||||
|
||||
void qsort_TYPE_noidx_big(TYPE* A, int32_t isize) {
|
||||
qsort( (void*) A, (size_t) isize, sizeof(TYPE), compare_TYPE_big);
|
||||
}
|
||||
#+END_SRC
|
||||
|
||||
* Fortran template
|
||||
|
||||
#+NAME:f_template
|
||||
#+BEGIN_SRC f90
|
||||
subroutine Lsort_big_c(A, iorder, isize) bind(C, name="qsort_TYPE_big")
|
||||
use iso_c_binding
|
||||
integer(c_int32_t), value :: isize
|
||||
integer(c_int32_t) :: iorder(isize)
|
||||
real (c_TYPE) :: A(isize)
|
||||
end subroutine Lsort_big_c
|
||||
|
||||
subroutine Lsort_noidx_big_c(A, isize) bind(C, name="qsort_TYPE_noidx_big")
|
||||
use iso_c_binding
|
||||
integer(c_int32_t), value :: isize
|
||||
real (c_TYPE) :: A(isize)
|
||||
end subroutine Lsort_noidx_big_c
|
||||
|
||||
#+END_SRC
|
||||
|
||||
#+NAME:f_template2
|
||||
#+BEGIN_SRC f90
|
||||
subroutine Lsort_big(A, iorder, isize)
|
||||
use qsort_module
|
||||
use iso_c_binding
|
||||
integer(c_int32_t) :: isize
|
||||
integer(c_int32_t) :: iorder(isize)
|
||||
real (c_TYPE) :: A(isize)
|
||||
call Lsort_big_c(A, iorder, isize)
|
||||
end subroutine Lsort_big
|
||||
|
||||
subroutine Lsort_noidx_big(A, isize)
|
||||
use iso_c_binding
|
||||
use qsort_module
|
||||
integer(c_int32_t) :: isize
|
||||
real (c_TYPE) :: A(isize)
|
||||
call Lsort_noidx_big_c(A, isize)
|
||||
end subroutine Lsort_noidx_big
|
||||
|
||||
#+END_SRC
|
||||
|
||||
* Python scripts for type replacements
|
||||
|
||||
#+NAME: replaced
|
||||
#+begin_src python :results output :noweb yes
|
||||
data = """
|
||||
<<c_template>>
|
||||
"""
|
||||
for typ in ["int16_t", "int32_t", "int64_t", "double", "float"]:
|
||||
print( data.replace("TYPE", typ).replace("_big", "") )
|
||||
print( data.replace("int32_t", "int64_t").replace("TYPE", typ) )
|
||||
#+end_src
|
||||
|
||||
#+NAME: replaced_f
|
||||
#+begin_src python :results output :noweb yes
|
||||
data = """
|
||||
<<f_template>>
|
||||
"""
|
||||
c1 = {
|
||||
"int16_t": "i2",
|
||||
"int32_t": "i",
|
||||
"int64_t": "i8",
|
||||
"double": "d",
|
||||
"float": ""
|
||||
}
|
||||
c2 = {
|
||||
"int16_t": "integer",
|
||||
"int32_t": "integer",
|
||||
"int64_t": "integer",
|
||||
"double": "real",
|
||||
"float": "real"
|
||||
}
|
||||
|
||||
for typ in ["int16_t", "int32_t", "int64_t", "double", "float"]:
|
||||
print( data.replace("real",c2[typ]).replace("L",c1[typ]).replace("TYPE", typ).replace("_big", "") )
|
||||
print( data.replace("real",c2[typ]).replace("L",c1[typ]).replace("int32_t", "int64_t").replace("TYPE", typ) )
|
||||
#+end_src
|
||||
|
||||
#+NAME: replaced_f2
|
||||
#+begin_src python :results output :noweb yes
|
||||
data = """
|
||||
<<f_template2>>
|
||||
"""
|
||||
c1 = {
|
||||
"int16_t": "i2",
|
||||
"int32_t": "i",
|
||||
"int64_t": "i8",
|
||||
"double": "d",
|
||||
"float": ""
|
||||
}
|
||||
c2 = {
|
||||
"int16_t": "integer",
|
||||
"int32_t": "integer",
|
||||
"int64_t": "integer",
|
||||
"double": "real",
|
||||
"float": "real"
|
||||
}
|
||||
|
||||
for typ in ["int16_t", "int32_t", "int64_t", "double", "float"]:
|
||||
print( data.replace("real",c2[typ]).replace("L",c1[typ]).replace("TYPE", typ).replace("_big", "") )
|
||||
print( data.replace("real",c2[typ]).replace("L",c1[typ]).replace("int32_t", "int64_t").replace("TYPE", typ) )
|
||||
#+end_src
|
||||
|
||||
* Generated C file
|
||||
|
||||
#+BEGIN_SRC c :comments link :tangle qsort.c :noweb yes
|
||||
#include <stdlib.h>
|
||||
#include <stdint.h>
|
||||
<<replaced()>>
|
||||
#+END_SRC
|
||||
|
||||
* Generated Fortran file
|
||||
|
||||
#+BEGIN_SRC f90 :tangle qsort_module.f90 :noweb yes
|
||||
module qsort_module
|
||||
use iso_c_binding
|
||||
|
||||
interface
|
||||
<<replaced_f()>>
|
||||
end interface
|
||||
|
||||
end module qsort_module
|
||||
|
||||
<<replaced_f2()>>
|
||||
|
||||
#+END_SRC
|
||||
|
||||
347
src/utils/qsort_module.f90
Normal file
347
src/utils/qsort_module.f90
Normal file
@ -0,0 +1,347 @@
|
||||
module qsort_module
|
||||
use iso_c_binding
|
||||
|
||||
interface
|
||||
|
||||
subroutine i2sort_c(A, iorder, isize) bind(C, name="qsort_int16_t")
|
||||
use iso_c_binding
|
||||
integer(c_int32_t), value :: isize
|
||||
integer(c_int32_t) :: iorder(isize)
|
||||
integer (c_int16_t) :: A(isize)
|
||||
end subroutine i2sort_c
|
||||
|
||||
subroutine i2sort_noidx_c(A, isize) bind(C, name="qsort_int16_t_noidx")
|
||||
use iso_c_binding
|
||||
integer(c_int32_t), value :: isize
|
||||
integer (c_int16_t) :: A(isize)
|
||||
end subroutine i2sort_noidx_c
|
||||
|
||||
|
||||
|
||||
subroutine i2sort_big_c(A, iorder, isize) bind(C, name="qsort_int16_t_big")
|
||||
use iso_c_binding
|
||||
integer(c_int64_t), value :: isize
|
||||
integer(c_int64_t) :: iorder(isize)
|
||||
integer (c_int16_t) :: A(isize)
|
||||
end subroutine i2sort_big_c
|
||||
|
||||
subroutine i2sort_noidx_big_c(A, isize) bind(C, name="qsort_int16_t_noidx_big")
|
||||
use iso_c_binding
|
||||
integer(c_int64_t), value :: isize
|
||||
integer (c_int16_t) :: A(isize)
|
||||
end subroutine i2sort_noidx_big_c
|
||||
|
||||
|
||||
|
||||
subroutine isort_c(A, iorder, isize) bind(C, name="qsort_int32_t")
|
||||
use iso_c_binding
|
||||
integer(c_int32_t), value :: isize
|
||||
integer(c_int32_t) :: iorder(isize)
|
||||
integer (c_int32_t) :: A(isize)
|
||||
end subroutine isort_c
|
||||
|
||||
subroutine isort_noidx_c(A, isize) bind(C, name="qsort_int32_t_noidx")
|
||||
use iso_c_binding
|
||||
integer(c_int32_t), value :: isize
|
||||
integer (c_int32_t) :: A(isize)
|
||||
end subroutine isort_noidx_c
|
||||
|
||||
|
||||
|
||||
subroutine isort_big_c(A, iorder, isize) bind(C, name="qsort_int32_t_big")
|
||||
use iso_c_binding
|
||||
integer(c_int64_t), value :: isize
|
||||
integer(c_int64_t) :: iorder(isize)
|
||||
integer (c_int32_t) :: A(isize)
|
||||
end subroutine isort_big_c
|
||||
|
||||
subroutine isort_noidx_big_c(A, isize) bind(C, name="qsort_int32_t_noidx_big")
|
||||
use iso_c_binding
|
||||
integer(c_int64_t), value :: isize
|
||||
integer (c_int32_t) :: A(isize)
|
||||
end subroutine isort_noidx_big_c
|
||||
|
||||
|
||||
|
||||
subroutine i8sort_c(A, iorder, isize) bind(C, name="qsort_int64_t")
|
||||
use iso_c_binding
|
||||
integer(c_int32_t), value :: isize
|
||||
integer(c_int32_t) :: iorder(isize)
|
||||
integer (c_int64_t) :: A(isize)
|
||||
end subroutine i8sort_c
|
||||
|
||||
subroutine i8sort_noidx_c(A, isize) bind(C, name="qsort_int64_t_noidx")
|
||||
use iso_c_binding
|
||||
integer(c_int32_t), value :: isize
|
||||
integer (c_int64_t) :: A(isize)
|
||||
end subroutine i8sort_noidx_c
|
||||
|
||||
|
||||
|
||||
subroutine i8sort_big_c(A, iorder, isize) bind(C, name="qsort_int64_t_big")
|
||||
use iso_c_binding
|
||||
integer(c_int64_t), value :: isize
|
||||
integer(c_int64_t) :: iorder(isize)
|
||||
integer (c_int64_t) :: A(isize)
|
||||
end subroutine i8sort_big_c
|
||||
|
||||
subroutine i8sort_noidx_big_c(A, isize) bind(C, name="qsort_int64_t_noidx_big")
|
||||
use iso_c_binding
|
||||
integer(c_int64_t), value :: isize
|
||||
integer (c_int64_t) :: A(isize)
|
||||
end subroutine i8sort_noidx_big_c
|
||||
|
||||
|
||||
|
||||
subroutine dsort_c(A, iorder, isize) bind(C, name="qsort_double")
|
||||
use iso_c_binding
|
||||
integer(c_int32_t), value :: isize
|
||||
integer(c_int32_t) :: iorder(isize)
|
||||
real (c_double) :: A(isize)
|
||||
end subroutine dsort_c
|
||||
|
||||
subroutine dsort_noidx_c(A, isize) bind(C, name="qsort_double_noidx")
|
||||
use iso_c_binding
|
||||
integer(c_int32_t), value :: isize
|
||||
real (c_double) :: A(isize)
|
||||
end subroutine dsort_noidx_c
|
||||
|
||||
|
||||
|
||||
subroutine dsort_big_c(A, iorder, isize) bind(C, name="qsort_double_big")
|
||||
use iso_c_binding
|
||||
integer(c_int64_t), value :: isize
|
||||
integer(c_int64_t) :: iorder(isize)
|
||||
real (c_double) :: A(isize)
|
||||
end subroutine dsort_big_c
|
||||
|
||||
subroutine dsort_noidx_big_c(A, isize) bind(C, name="qsort_double_noidx_big")
|
||||
use iso_c_binding
|
||||
integer(c_int64_t), value :: isize
|
||||
real (c_double) :: A(isize)
|
||||
end subroutine dsort_noidx_big_c
|
||||
|
||||
|
||||
|
||||
subroutine sort_c(A, iorder, isize) bind(C, name="qsort_float")
|
||||
use iso_c_binding
|
||||
integer(c_int32_t), value :: isize
|
||||
integer(c_int32_t) :: iorder(isize)
|
||||
real (c_float) :: A(isize)
|
||||
end subroutine sort_c
|
||||
|
||||
subroutine sort_noidx_c(A, isize) bind(C, name="qsort_float_noidx")
|
||||
use iso_c_binding
|
||||
integer(c_int32_t), value :: isize
|
||||
real (c_float) :: A(isize)
|
||||
end subroutine sort_noidx_c
|
||||
|
||||
|
||||
|
||||
subroutine sort_big_c(A, iorder, isize) bind(C, name="qsort_float_big")
|
||||
use iso_c_binding
|
||||
integer(c_int64_t), value :: isize
|
||||
integer(c_int64_t) :: iorder(isize)
|
||||
real (c_float) :: A(isize)
|
||||
end subroutine sort_big_c
|
||||
|
||||
subroutine sort_noidx_big_c(A, isize) bind(C, name="qsort_float_noidx_big")
|
||||
use iso_c_binding
|
||||
integer(c_int64_t), value :: isize
|
||||
real (c_float) :: A(isize)
|
||||
end subroutine sort_noidx_big_c
|
||||
|
||||
|
||||
|
||||
end interface
|
||||
|
||||
end module qsort_module
|
||||
|
||||
|
||||
subroutine i2sort(A, iorder, isize)
|
||||
use qsort_module
|
||||
use iso_c_binding
|
||||
integer(c_int32_t) :: isize
|
||||
integer(c_int32_t) :: iorder(isize)
|
||||
integer (c_int16_t) :: A(isize)
|
||||
call i2sort_c(A, iorder, isize)
|
||||
end subroutine i2sort
|
||||
|
||||
subroutine i2sort_noidx(A, isize)
|
||||
use iso_c_binding
|
||||
use qsort_module
|
||||
integer(c_int32_t) :: isize
|
||||
integer (c_int16_t) :: A(isize)
|
||||
call i2sort_noidx_c(A, isize)
|
||||
end subroutine i2sort_noidx
|
||||
|
||||
|
||||
|
||||
subroutine i2sort_big(A, iorder, isize)
|
||||
use qsort_module
|
||||
use iso_c_binding
|
||||
integer(c_int64_t) :: isize
|
||||
integer(c_int64_t) :: iorder(isize)
|
||||
integer (c_int16_t) :: A(isize)
|
||||
call i2sort_big_c(A, iorder, isize)
|
||||
end subroutine i2sort_big
|
||||
|
||||
subroutine i2sort_noidx_big(A, isize)
|
||||
use iso_c_binding
|
||||
use qsort_module
|
||||
integer(c_int64_t) :: isize
|
||||
integer (c_int16_t) :: A(isize)
|
||||
call i2sort_noidx_big_c(A, isize)
|
||||
end subroutine i2sort_noidx_big
|
||||
|
||||
|
||||
|
||||
subroutine isort(A, iorder, isize)
|
||||
use qsort_module
|
||||
use iso_c_binding
|
||||
integer(c_int32_t) :: isize
|
||||
integer(c_int32_t) :: iorder(isize)
|
||||
integer (c_int32_t) :: A(isize)
|
||||
call isort_c(A, iorder, isize)
|
||||
end subroutine isort
|
||||
|
||||
subroutine isort_noidx(A, isize)
|
||||
use iso_c_binding
|
||||
use qsort_module
|
||||
integer(c_int32_t) :: isize
|
||||
integer (c_int32_t) :: A(isize)
|
||||
call isort_noidx_c(A, isize)
|
||||
end subroutine isort_noidx
|
||||
|
||||
|
||||
|
||||
subroutine isort_big(A, iorder, isize)
|
||||
use qsort_module
|
||||
use iso_c_binding
|
||||
integer(c_int64_t) :: isize
|
||||
integer(c_int64_t) :: iorder(isize)
|
||||
integer (c_int32_t) :: A(isize)
|
||||
call isort_big_c(A, iorder, isize)
|
||||
end subroutine isort_big
|
||||
|
||||
subroutine isort_noidx_big(A, isize)
|
||||
use iso_c_binding
|
||||
use qsort_module
|
||||
integer(c_int64_t) :: isize
|
||||
integer (c_int32_t) :: A(isize)
|
||||
call isort_noidx_big_c(A, isize)
|
||||
end subroutine isort_noidx_big
|
||||
|
||||
|
||||
|
||||
subroutine i8sort(A, iorder, isize)
|
||||
use qsort_module
|
||||
use iso_c_binding
|
||||
integer(c_int32_t) :: isize
|
||||
integer(c_int32_t) :: iorder(isize)
|
||||
integer (c_int64_t) :: A(isize)
|
||||
call i8sort_c(A, iorder, isize)
|
||||
end subroutine i8sort
|
||||
|
||||
subroutine i8sort_noidx(A, isize)
|
||||
use iso_c_binding
|
||||
use qsort_module
|
||||
integer(c_int32_t) :: isize
|
||||
integer (c_int64_t) :: A(isize)
|
||||
call i8sort_noidx_c(A, isize)
|
||||
end subroutine i8sort_noidx
|
||||
|
||||
|
||||
|
||||
subroutine i8sort_big(A, iorder, isize)
|
||||
use qsort_module
|
||||
use iso_c_binding
|
||||
integer(c_int64_t) :: isize
|
||||
integer(c_int64_t) :: iorder(isize)
|
||||
integer (c_int64_t) :: A(isize)
|
||||
call i8sort_big_c(A, iorder, isize)
|
||||
end subroutine i8sort_big
|
||||
|
||||
subroutine i8sort_noidx_big(A, isize)
|
||||
use iso_c_binding
|
||||
use qsort_module
|
||||
integer(c_int64_t) :: isize
|
||||
integer (c_int64_t) :: A(isize)
|
||||
call i8sort_noidx_big_c(A, isize)
|
||||
end subroutine i8sort_noidx_big
|
||||
|
||||
|
||||
|
||||
subroutine dsort(A, iorder, isize)
|
||||
use qsort_module
|
||||
use iso_c_binding
|
||||
integer(c_int32_t) :: isize
|
||||
integer(c_int32_t) :: iorder(isize)
|
||||
real (c_double) :: A(isize)
|
||||
call dsort_c(A, iorder, isize)
|
||||
end subroutine dsort
|
||||
|
||||
subroutine dsort_noidx(A, isize)
|
||||
use iso_c_binding
|
||||
use qsort_module
|
||||
integer(c_int32_t) :: isize
|
||||
real (c_double) :: A(isize)
|
||||
call dsort_noidx_c(A, isize)
|
||||
end subroutine dsort_noidx
|
||||
|
||||
|
||||
|
||||
subroutine dsort_big(A, iorder, isize)
|
||||
use qsort_module
|
||||
use iso_c_binding
|
||||
integer(c_int64_t) :: isize
|
||||
integer(c_int64_t) :: iorder(isize)
|
||||
real (c_double) :: A(isize)
|
||||
call dsort_big_c(A, iorder, isize)
|
||||
end subroutine dsort_big
|
||||
|
||||
subroutine dsort_noidx_big(A, isize)
|
||||
use iso_c_binding
|
||||
use qsort_module
|
||||
integer(c_int64_t) :: isize
|
||||
real (c_double) :: A(isize)
|
||||
call dsort_noidx_big_c(A, isize)
|
||||
end subroutine dsort_noidx_big
|
||||
|
||||
|
||||
|
||||
subroutine sort(A, iorder, isize)
|
||||
use qsort_module
|
||||
use iso_c_binding
|
||||
integer(c_int32_t) :: isize
|
||||
integer(c_int32_t) :: iorder(isize)
|
||||
real (c_float) :: A(isize)
|
||||
call sort_c(A, iorder, isize)
|
||||
end subroutine sort
|
||||
|
||||
subroutine sort_noidx(A, isize)
|
||||
use iso_c_binding
|
||||
use qsort_module
|
||||
integer(c_int32_t) :: isize
|
||||
real (c_float) :: A(isize)
|
||||
call sort_noidx_c(A, isize)
|
||||
end subroutine sort_noidx
|
||||
|
||||
|
||||
|
||||
subroutine sort_big(A, iorder, isize)
|
||||
use qsort_module
|
||||
use iso_c_binding
|
||||
integer(c_int64_t) :: isize
|
||||
integer(c_int64_t) :: iorder(isize)
|
||||
real (c_float) :: A(isize)
|
||||
call sort_big_c(A, iorder, isize)
|
||||
end subroutine sort_big
|
||||
|
||||
subroutine sort_noidx_big(A, isize)
|
||||
use iso_c_binding
|
||||
use qsort_module
|
||||
integer(c_int64_t) :: isize
|
||||
real (c_float) :: A(isize)
|
||||
call sort_noidx_big_c(A, isize)
|
||||
end subroutine sort_noidx_big
|
||||
@ -1,222 +1,4 @@
|
||||
BEGIN_TEMPLATE
|
||||
subroutine insertion_$Xsort (x,iorder,isize)
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Sort array x(isize) using the insertion sort algorithm.
|
||||
! iorder in input should be (1,2,3,...,isize), and in output
|
||||
! contains the new order of the elements.
|
||||
END_DOC
|
||||
integer,intent(in) :: isize
|
||||
$type,intent(inout) :: x(isize)
|
||||
integer,intent(inout) :: iorder(isize)
|
||||
$type :: xtmp
|
||||
integer :: i, i0, j, jmax
|
||||
|
||||
do i=2,isize
|
||||
xtmp = x(i)
|
||||
i0 = iorder(i)
|
||||
j=i-1
|
||||
do while (j>0)
|
||||
if ((x(j) <= xtmp)) exit
|
||||
x(j+1) = x(j)
|
||||
iorder(j+1) = iorder(j)
|
||||
j=j-1
|
||||
enddo
|
||||
x(j+1) = xtmp
|
||||
iorder(j+1) = i0
|
||||
enddo
|
||||
end subroutine insertion_$Xsort
|
||||
|
||||
subroutine quick_$Xsort(x, iorder, isize)
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Sort array x(isize) using the quicksort algorithm.
|
||||
! iorder in input should be (1,2,3,...,isize), and in output
|
||||
! contains the new order of the elements.
|
||||
END_DOC
|
||||
integer,intent(in) :: isize
|
||||
$type,intent(inout) :: x(isize)
|
||||
integer,intent(inout) :: iorder(isize)
|
||||
integer, external :: omp_get_num_threads
|
||||
call rec_$X_quicksort(x,iorder,isize,1,isize,nproc)
|
||||
end
|
||||
|
||||
recursive subroutine rec_$X_quicksort(x, iorder, isize, first, last, level)
|
||||
implicit none
|
||||
integer, intent(in) :: isize, first, last, level
|
||||
integer,intent(inout) :: iorder(isize)
|
||||
$type, intent(inout) :: x(isize)
|
||||
$type :: c, tmp
|
||||
integer :: itmp
|
||||
integer :: i, j
|
||||
|
||||
if(isize<2)return
|
||||
|
||||
c = x( shiftr(first+last,1) )
|
||||
i = first
|
||||
j = last
|
||||
do
|
||||
do while (x(i) < c)
|
||||
i=i+1
|
||||
end do
|
||||
do while (c < x(j))
|
||||
j=j-1
|
||||
end do
|
||||
if (i >= j) exit
|
||||
tmp = x(i)
|
||||
x(i) = x(j)
|
||||
x(j) = tmp
|
||||
itmp = iorder(i)
|
||||
iorder(i) = iorder(j)
|
||||
iorder(j) = itmp
|
||||
i=i+1
|
||||
j=j-1
|
||||
enddo
|
||||
if ( ((i-first <= 10000).and.(last-j <= 10000)).or.(level<=0) ) then
|
||||
if (first < i-1) then
|
||||
call rec_$X_quicksort(x, iorder, isize, first, i-1,level/2)
|
||||
endif
|
||||
if (j+1 < last) then
|
||||
call rec_$X_quicksort(x, iorder, isize, j+1, last,level/2)
|
||||
endif
|
||||
else
|
||||
if (first < i-1) then
|
||||
call rec_$X_quicksort(x, iorder, isize, first, i-1,level/2)
|
||||
endif
|
||||
if (j+1 < last) then
|
||||
call rec_$X_quicksort(x, iorder, isize, j+1, last,level/2)
|
||||
endif
|
||||
endif
|
||||
end
|
||||
|
||||
subroutine heap_$Xsort(x,iorder,isize)
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Sort array x(isize) using the heap sort algorithm.
|
||||
! iorder in input should be (1,2,3,...,isize), and in output
|
||||
! contains the new order of the elements.
|
||||
END_DOC
|
||||
integer,intent(in) :: isize
|
||||
$type,intent(inout) :: x(isize)
|
||||
integer,intent(inout) :: iorder(isize)
|
||||
|
||||
integer :: i, k, j, l, i0
|
||||
$type :: xtemp
|
||||
|
||||
l = isize/2+1
|
||||
k = isize
|
||||
do while (.True.)
|
||||
if (l>1) then
|
||||
l=l-1
|
||||
xtemp = x(l)
|
||||
i0 = iorder(l)
|
||||
else
|
||||
xtemp = x(k)
|
||||
i0 = iorder(k)
|
||||
x(k) = x(1)
|
||||
iorder(k) = iorder(1)
|
||||
k = k-1
|
||||
if (k == 1) then
|
||||
x(1) = xtemp
|
||||
iorder(1) = i0
|
||||
exit
|
||||
endif
|
||||
endif
|
||||
i=l
|
||||
j = shiftl(l,1)
|
||||
do while (j<k)
|
||||
if ( x(j) < x(j+1) ) then
|
||||
j=j+1
|
||||
endif
|
||||
if (xtemp < x(j)) then
|
||||
x(i) = x(j)
|
||||
iorder(i) = iorder(j)
|
||||
i = j
|
||||
j = shiftl(j,1)
|
||||
else
|
||||
j = k+1
|
||||
endif
|
||||
enddo
|
||||
if (j==k) then
|
||||
if (xtemp < x(j)) then
|
||||
x(i) = x(j)
|
||||
iorder(i) = iorder(j)
|
||||
i = j
|
||||
j = shiftl(j,1)
|
||||
else
|
||||
j = k+1
|
||||
endif
|
||||
endif
|
||||
x(i) = xtemp
|
||||
iorder(i) = i0
|
||||
enddo
|
||||
end subroutine heap_$Xsort
|
||||
|
||||
subroutine heap_$Xsort_big(x,iorder,isize)
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Sort array x(isize) using the heap sort algorithm.
|
||||
! iorder in input should be (1,2,3,...,isize), and in output
|
||||
! contains the new order of the elements.
|
||||
! This is a version for very large arrays where the indices need
|
||||
! to be in integer*8 format
|
||||
END_DOC
|
||||
integer*8,intent(in) :: isize
|
||||
$type,intent(inout) :: x(isize)
|
||||
integer*8,intent(inout) :: iorder(isize)
|
||||
|
||||
integer*8 :: i, k, j, l, i0
|
||||
$type :: xtemp
|
||||
|
||||
l = isize/2+1
|
||||
k = isize
|
||||
do while (.True.)
|
||||
if (l>1) then
|
||||
l=l-1
|
||||
xtemp = x(l)
|
||||
i0 = iorder(l)
|
||||
else
|
||||
xtemp = x(k)
|
||||
i0 = iorder(k)
|
||||
x(k) = x(1)
|
||||
iorder(k) = iorder(1)
|
||||
k = k-1
|
||||
if (k == 1) then
|
||||
x(1) = xtemp
|
||||
iorder(1) = i0
|
||||
exit
|
||||
endif
|
||||
endif
|
||||
i=l
|
||||
j = shiftl(l,1)
|
||||
do while (j<k)
|
||||
if ( x(j) < x(j+1) ) then
|
||||
j=j+1
|
||||
endif
|
||||
if (xtemp < x(j)) then
|
||||
x(i) = x(j)
|
||||
iorder(i) = iorder(j)
|
||||
i = j
|
||||
j = shiftl(j,1)
|
||||
else
|
||||
j = k+1
|
||||
endif
|
||||
enddo
|
||||
if (j==k) then
|
||||
if (xtemp < x(j)) then
|
||||
x(i) = x(j)
|
||||
iorder(i) = iorder(j)
|
||||
i = j
|
||||
j = shiftl(j,1)
|
||||
else
|
||||
j = k+1
|
||||
endif
|
||||
endif
|
||||
x(i) = xtemp
|
||||
iorder(i) = i0
|
||||
enddo
|
||||
|
||||
end subroutine heap_$Xsort_big
|
||||
|
||||
subroutine sorted_$Xnumber(x,isize,n)
|
||||
implicit none
|
||||
@ -250,222 +32,6 @@ SUBST [ X, type ]
|
||||
END_TEMPLATE
|
||||
|
||||
|
||||
!---------------------- INTEL
|
||||
IRP_IF INTEL
|
||||
|
||||
BEGIN_TEMPLATE
|
||||
subroutine $Xsort(x,iorder,isize)
|
||||
use intel
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Sort array x(isize).
|
||||
! iorder in input should be (1,2,3,...,isize), and in output
|
||||
! contains the new order of the elements.
|
||||
END_DOC
|
||||
integer,intent(in) :: isize
|
||||
$type,intent(inout) :: x(isize)
|
||||
integer,intent(inout) :: iorder(isize)
|
||||
integer :: n
|
||||
character, allocatable :: tmp(:)
|
||||
if (isize < 2) return
|
||||
call ippsSortRadixIndexGetBufferSize(isize, $ippsz, n)
|
||||
allocate(tmp(n))
|
||||
call ippsSortRadixIndexAscend_$ityp(x, $n, iorder, isize, tmp)
|
||||
deallocate(tmp)
|
||||
iorder(1:isize) = iorder(1:isize)+1
|
||||
call $Xset_order(x,iorder,isize)
|
||||
end
|
||||
|
||||
subroutine $Xsort_noidx(x,isize)
|
||||
use intel
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Sort array x(isize).
|
||||
! iorder in input should be (1,2,3,...,isize), and in output
|
||||
! contains the new order of the elements.
|
||||
END_DOC
|
||||
integer,intent(in) :: isize
|
||||
$type,intent(inout) :: x(isize)
|
||||
integer :: n
|
||||
character, allocatable :: tmp(:)
|
||||
if (isize < 2) return
|
||||
call ippsSortRadixIndexGetBufferSize(isize, $ippsz, n)
|
||||
allocate(tmp(n))
|
||||
call ippsSortRadixAscend_$ityp_I(x, isize, tmp)
|
||||
deallocate(tmp)
|
||||
end
|
||||
|
||||
SUBST [ X, type, ityp, n, ippsz ]
|
||||
; real ; 32f ; 4 ; 13 ;;
|
||||
i ; integer ; 32s ; 4 ; 11 ;;
|
||||
i2 ; integer*2 ; 16s ; 2 ; 7 ;;
|
||||
END_TEMPLATE
|
||||
|
||||
BEGIN_TEMPLATE
|
||||
|
||||
subroutine $Xsort(x,iorder,isize)
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Sort array x(isize).
|
||||
! iorder in input should be (1,2,3,...,isize), and in output
|
||||
! contains the new order of the elements.
|
||||
END_DOC
|
||||
integer,intent(in) :: isize
|
||||
$type,intent(inout) :: x(isize)
|
||||
integer,intent(inout) :: iorder(isize)
|
||||
integer :: n
|
||||
if (isize < 2) then
|
||||
return
|
||||
endif
|
||||
! call sorted_$Xnumber(x,isize,n)
|
||||
! if (isize == n) then
|
||||
! return
|
||||
! endif
|
||||
if ( isize < 32) then
|
||||
call insertion_$Xsort(x,iorder,isize)
|
||||
else
|
||||
! call heap_$Xsort(x,iorder,isize)
|
||||
call quick_$Xsort(x,iorder,isize)
|
||||
endif
|
||||
end subroutine $Xsort
|
||||
|
||||
SUBST [ X, type ]
|
||||
d ; double precision ;;
|
||||
END_TEMPLATE
|
||||
|
||||
BEGIN_TEMPLATE
|
||||
|
||||
subroutine $Xsort(x,iorder,isize)
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Sort array x(isize).
|
||||
! iorder in input should be (1,2,3,...,isize), and in output
|
||||
! contains the new order of the elements.
|
||||
END_DOC
|
||||
integer,intent(in) :: isize
|
||||
$type,intent(inout) :: x(isize)
|
||||
integer,intent(inout) :: iorder(isize)
|
||||
integer :: n
|
||||
if (isize < 2) then
|
||||
return
|
||||
endif
|
||||
call sorted_$Xnumber(x,isize,n)
|
||||
if (isize == n) then
|
||||
return
|
||||
endif
|
||||
if ( isize < 32) then
|
||||
call insertion_$Xsort(x,iorder,isize)
|
||||
else
|
||||
! call $Xradix_sort(x,iorder,isize,-1)
|
||||
call quick_$Xsort(x,iorder,isize)
|
||||
endif
|
||||
end subroutine $Xsort
|
||||
|
||||
SUBST [ X, type ]
|
||||
i8 ; integer*8 ;;
|
||||
END_TEMPLATE
|
||||
|
||||
!---------------------- END INTEL
|
||||
IRP_ELSE
|
||||
!---------------------- NON-INTEL
|
||||
BEGIN_TEMPLATE
|
||||
|
||||
subroutine $Xsort_noidx(x,isize)
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Sort array x(isize).
|
||||
END_DOC
|
||||
integer,intent(in) :: isize
|
||||
$type,intent(inout) :: x(isize)
|
||||
integer, allocatable :: iorder(:)
|
||||
integer :: i
|
||||
allocate(iorder(isize))
|
||||
do i=1,isize
|
||||
iorder(i)=i
|
||||
enddo
|
||||
call $Xsort(x,iorder,isize)
|
||||
deallocate(iorder)
|
||||
end subroutine $Xsort_noidx
|
||||
|
||||
SUBST [ X, type ]
|
||||
; real ;;
|
||||
d ; double precision ;;
|
||||
i ; integer ;;
|
||||
i8 ; integer*8 ;;
|
||||
i2 ; integer*2 ;;
|
||||
END_TEMPLATE
|
||||
|
||||
BEGIN_TEMPLATE
|
||||
|
||||
subroutine $Xsort(x,iorder,isize)
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Sort array x(isize).
|
||||
! iorder in input should be (1,2,3,...,isize), and in output
|
||||
! contains the new order of the elements.
|
||||
END_DOC
|
||||
integer,intent(in) :: isize
|
||||
$type,intent(inout) :: x(isize)
|
||||
integer,intent(inout) :: iorder(isize)
|
||||
integer :: n
|
||||
if (isize < 2) then
|
||||
return
|
||||
endif
|
||||
! call sorted_$Xnumber(x,isize,n)
|
||||
! if (isize == n) then
|
||||
! return
|
||||
! endif
|
||||
if ( isize < 32) then
|
||||
call insertion_$Xsort(x,iorder,isize)
|
||||
else
|
||||
! call heap_$Xsort(x,iorder,isize)
|
||||
call quick_$Xsort(x,iorder,isize)
|
||||
endif
|
||||
end subroutine $Xsort
|
||||
|
||||
SUBST [ X, type ]
|
||||
; real ;;
|
||||
d ; double precision ;;
|
||||
END_TEMPLATE
|
||||
|
||||
BEGIN_TEMPLATE
|
||||
|
||||
subroutine $Xsort(x,iorder,isize)
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Sort array x(isize).
|
||||
! iorder in input should be (1,2,3,...,isize), and in output
|
||||
! contains the new order of the elements.
|
||||
END_DOC
|
||||
integer,intent(in) :: isize
|
||||
$type,intent(inout) :: x(isize)
|
||||
integer,intent(inout) :: iorder(isize)
|
||||
integer :: n
|
||||
if (isize < 2) then
|
||||
return
|
||||
endif
|
||||
call sorted_$Xnumber(x,isize,n)
|
||||
if (isize == n) then
|
||||
return
|
||||
endif
|
||||
if ( isize < 32) then
|
||||
call insertion_$Xsort(x,iorder,isize)
|
||||
else
|
||||
! call $Xradix_sort(x,iorder,isize,-1)
|
||||
call quick_$Xsort(x,iorder,isize)
|
||||
endif
|
||||
end subroutine $Xsort
|
||||
|
||||
SUBST [ X, type ]
|
||||
i ; integer ;;
|
||||
i8 ; integer*8 ;;
|
||||
i2 ; integer*2 ;;
|
||||
END_TEMPLATE
|
||||
|
||||
IRP_ENDIF
|
||||
!---------------------- END NON-INTEL
|
||||
|
||||
|
||||
|
||||
BEGIN_TEMPLATE
|
||||
subroutine $Xset_order(x,iorder,isize)
|
||||
@ -491,47 +57,6 @@ BEGIN_TEMPLATE
|
||||
deallocate(xtmp)
|
||||
end
|
||||
|
||||
SUBST [ X, type ]
|
||||
; real ;;
|
||||
d ; double precision ;;
|
||||
i ; integer ;;
|
||||
i8; integer*8 ;;
|
||||
i2; integer*2 ;;
|
||||
END_TEMPLATE
|
||||
|
||||
|
||||
BEGIN_TEMPLATE
|
||||
subroutine insertion_$Xsort_big (x,iorder,isize)
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Sort array x(isize) using the insertion sort algorithm.
|
||||
! iorder in input should be (1,2,3,...,isize), and in output
|
||||
! contains the new order of the elements.
|
||||
! This is a version for very large arrays where the indices need
|
||||
! to be in integer*8 format
|
||||
END_DOC
|
||||
integer*8,intent(in) :: isize
|
||||
$type,intent(inout) :: x(isize)
|
||||
integer*8,intent(inout) :: iorder(isize)
|
||||
$type :: xtmp
|
||||
integer*8 :: i, i0, j, jmax
|
||||
|
||||
do i=2_8,isize
|
||||
xtmp = x(i)
|
||||
i0 = iorder(i)
|
||||
j = i-1_8
|
||||
do while (j>0_8)
|
||||
if (x(j)<=xtmp) exit
|
||||
x(j+1_8) = x(j)
|
||||
iorder(j+1_8) = iorder(j)
|
||||
j = j-1_8
|
||||
enddo
|
||||
x(j+1_8) = xtmp
|
||||
iorder(j+1_8) = i0
|
||||
enddo
|
||||
|
||||
end subroutine insertion_$Xsort_big
|
||||
|
||||
subroutine $Xset_order_big(x,iorder,isize)
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
@ -565,223 +90,3 @@ SUBST [ X, type ]
|
||||
END_TEMPLATE
|
||||
|
||||
|
||||
BEGIN_TEMPLATE
|
||||
|
||||
recursive subroutine $Xradix_sort$big(x,iorder,isize,iradix)
|
||||
implicit none
|
||||
|
||||
BEGIN_DOC
|
||||
! Sort integer array x(isize) using the radix sort algorithm.
|
||||
! iorder in input should be (1,2,3,...,isize), and in output
|
||||
! contains the new order of the elements.
|
||||
! iradix should be -1 in input.
|
||||
END_DOC
|
||||
integer*$int_type, intent(in) :: isize
|
||||
integer*$int_type, intent(inout) :: iorder(isize)
|
||||
integer*$type, intent(inout) :: x(isize)
|
||||
integer, intent(in) :: iradix
|
||||
integer :: iradix_new
|
||||
integer*$type, allocatable :: x2(:), x1(:)
|
||||
integer*$type :: i4 ! data type
|
||||
integer*$int_type, allocatable :: iorder1(:),iorder2(:)
|
||||
integer*$int_type :: i0, i1, i2, i3, i ! index type
|
||||
integer*$type :: mask
|
||||
integer :: err
|
||||
!DIR$ ATTRIBUTES ALIGN : 128 :: iorder1,iorder2, x2, x1
|
||||
|
||||
if (isize < 2) then
|
||||
return
|
||||
endif
|
||||
|
||||
if (iradix == -1) then ! Sort Positive and negative
|
||||
|
||||
allocate(x1(isize),iorder1(isize), x2(isize),iorder2(isize),stat=err)
|
||||
if (err /= 0) then
|
||||
print *, irp_here, ': Unable to allocate arrays'
|
||||
stop
|
||||
endif
|
||||
|
||||
i1=1_$int_type
|
||||
i2=1_$int_type
|
||||
do i=1_$int_type,isize
|
||||
if (x(i) < 0_$type) then
|
||||
iorder1(i1) = iorder(i)
|
||||
x1(i1) = -x(i)
|
||||
i1 = i1+1_$int_type
|
||||
else
|
||||
iorder2(i2) = iorder(i)
|
||||
x2(i2) = x(i)
|
||||
i2 = i2+1_$int_type
|
||||
endif
|
||||
enddo
|
||||
i1=i1-1_$int_type
|
||||
i2=i2-1_$int_type
|
||||
|
||||
do i=1_$int_type,i2
|
||||
iorder(i1+i) = iorder2(i)
|
||||
x(i1+i) = x2(i)
|
||||
enddo
|
||||
deallocate(x2,iorder2,stat=err)
|
||||
if (err /= 0) then
|
||||
print *, irp_here, ': Unable to deallocate arrays x2, iorder2'
|
||||
stop
|
||||
endif
|
||||
|
||||
|
||||
if (i1 > 1_$int_type) then
|
||||
call $Xradix_sort$big(x1,iorder1,i1,-2)
|
||||
do i=1_$int_type,i1
|
||||
x(i) = -x1(1_$int_type+i1-i)
|
||||
iorder(i) = iorder1(1_$int_type+i1-i)
|
||||
enddo
|
||||
endif
|
||||
|
||||
if (i2>1_$int_type) then
|
||||
call $Xradix_sort$big(x(i1+1_$int_type),iorder(i1+1_$int_type),i2,-2)
|
||||
endif
|
||||
|
||||
deallocate(x1,iorder1,stat=err)
|
||||
if (err /= 0) then
|
||||
print *, irp_here, ': Unable to deallocate arrays x1, iorder1'
|
||||
stop
|
||||
endif
|
||||
return
|
||||
|
||||
else if (iradix == -2) then ! Positive
|
||||
|
||||
! Find most significant bit
|
||||
|
||||
i0 = 0_$int_type
|
||||
i4 = maxval(x)
|
||||
|
||||
iradix_new = max($integer_size-1-leadz(i4),1)
|
||||
mask = ibset(0_$type,iradix_new)
|
||||
|
||||
allocate(x1(isize),iorder1(isize), x2(isize),iorder2(isize),stat=err)
|
||||
if (err /= 0) then
|
||||
print *, irp_here, ': Unable to allocate arrays'
|
||||
stop
|
||||
endif
|
||||
|
||||
i1=1_$int_type
|
||||
i2=1_$int_type
|
||||
|
||||
do i=1_$int_type,isize
|
||||
if (iand(mask,x(i)) == 0_$type) then
|
||||
iorder1(i1) = iorder(i)
|
||||
x1(i1) = x(i)
|
||||
i1 = i1+1_$int_type
|
||||
else
|
||||
iorder2(i2) = iorder(i)
|
||||
x2(i2) = x(i)
|
||||
i2 = i2+1_$int_type
|
||||
endif
|
||||
enddo
|
||||
i1=i1-1_$int_type
|
||||
i2=i2-1_$int_type
|
||||
|
||||
do i=1_$int_type,i1
|
||||
iorder(i0+i) = iorder1(i)
|
||||
x(i0+i) = x1(i)
|
||||
enddo
|
||||
i0 = i0+i1
|
||||
i3 = i0
|
||||
deallocate(x1,iorder1,stat=err)
|
||||
if (err /= 0) then
|
||||
print *, irp_here, ': Unable to deallocate arrays x1, iorder1'
|
||||
stop
|
||||
endif
|
||||
|
||||
|
||||
do i=1_$int_type,i2
|
||||
iorder(i0+i) = iorder2(i)
|
||||
x(i0+i) = x2(i)
|
||||
enddo
|
||||
i0 = i0+i2
|
||||
deallocate(x2,iorder2,stat=err)
|
||||
if (err /= 0) then
|
||||
print *, irp_here, ': Unable to deallocate arrays x2, iorder2'
|
||||
stop
|
||||
endif
|
||||
|
||||
|
||||
if (i3>1_$int_type) then
|
||||
call $Xradix_sort$big(x,iorder,i3,iradix_new-1)
|
||||
endif
|
||||
|
||||
if (isize-i3>1_$int_type) then
|
||||
call $Xradix_sort$big(x(i3+1_$int_type),iorder(i3+1_$int_type),isize-i3,iradix_new-1)
|
||||
endif
|
||||
|
||||
return
|
||||
endif
|
||||
|
||||
ASSERT (iradix >= 0)
|
||||
|
||||
if (isize < 48) then
|
||||
call insertion_$Xsort$big(x,iorder,isize)
|
||||
return
|
||||
endif
|
||||
|
||||
|
||||
allocate(x2(isize),iorder2(isize),stat=err)
|
||||
if (err /= 0) then
|
||||
print *, irp_here, ': Unable to allocate arrays x1, iorder1'
|
||||
stop
|
||||
endif
|
||||
|
||||
|
||||
mask = ibset(0_$type,iradix)
|
||||
i0=1_$int_type
|
||||
i1=1_$int_type
|
||||
|
||||
do i=1_$int_type,isize
|
||||
if (iand(mask,x(i)) == 0_$type) then
|
||||
iorder(i0) = iorder(i)
|
||||
x(i0) = x(i)
|
||||
i0 = i0+1_$int_type
|
||||
else
|
||||
iorder2(i1) = iorder(i)
|
||||
x2(i1) = x(i)
|
||||
i1 = i1+1_$int_type
|
||||
endif
|
||||
enddo
|
||||
i0=i0-1_$int_type
|
||||
i1=i1-1_$int_type
|
||||
|
||||
do i=1_$int_type,i1
|
||||
iorder(i0+i) = iorder2(i)
|
||||
x(i0+i) = x2(i)
|
||||
enddo
|
||||
|
||||
deallocate(x2,iorder2,stat=err)
|
||||
if (err /= 0) then
|
||||
print *, irp_here, ': Unable to allocate arrays x2, iorder2'
|
||||
stop
|
||||
endif
|
||||
|
||||
|
||||
if (iradix == 0) then
|
||||
return
|
||||
endif
|
||||
|
||||
|
||||
if (i1>1_$int_type) then
|
||||
call $Xradix_sort$big(x(i0+1_$int_type),iorder(i0+1_$int_type),i1,iradix-1)
|
||||
endif
|
||||
if (i0>1) then
|
||||
call $Xradix_sort$big(x,iorder,i0,iradix-1)
|
||||
endif
|
||||
|
||||
end
|
||||
|
||||
SUBST [ X, type, integer_size, is_big, big, int_type ]
|
||||
i ; 4 ; 32 ; .False. ; ; 4 ;;
|
||||
i8 ; 8 ; 64 ; .False. ; ; 4 ;;
|
||||
i2 ; 2 ; 16 ; .False. ; ; 4 ;;
|
||||
i ; 4 ; 32 ; .True. ; _big ; 8 ;;
|
||||
i8 ; 8 ; 64 ; .True. ; _big ; 8 ;;
|
||||
END_TEMPLATE
|
||||
|
||||
|
||||
|
||||
|
||||
Loading…
Reference in New Issue
Block a user