Skip to content
Snippets Groups Projects
Commit 1aaea3d6 authored by Giovanni La Mura's avatar Giovanni La Mura
Browse files

Merge branch 'link_with_lapacke' into 'master' 

BUG FIX: reshape C1 common rc field to configuration number

See merge request giacomo.mulas/np_tmcode!23
parents 4d4a311d a5ff64a7
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
src/libnptm/Commons.cpp
+ 208
− 194
View file @ 1aaea3d6
......@@ -19,234 +19,248 @@
using namespace std;
C1::C1(int ns, int l_max, int *_nshl, int *_iog) {
nlmmt = 2 * (l_max * (l_max + 2));
nsph = ns;
lm = l_max;
rmi = new complex<double>*[lm];
rei = new complex<double>*[lm];
for (int ri = 0; ri < lm; ri++) {
rmi[ri] = new complex<double>[nsph]();
rei[ri] = new complex<double>[nsph]();
}
w = new complex<double>*[nlmmt];
for (int wi = 0; wi < nlmmt; wi++) w[wi] = new complex<double>[4]();
vints = new complex<double>*[nsph];
rc = new double*[nsph];
nshl = new int[nsph]();
iog = new int[nsph]();
for (int vi = 0; vi < nsph; vi++) {
rc[vi] = new double[_nshl[vi]]();
vints[vi] = new complex<double>[16]();
nshl[vi] = _nshl[vi];
iog[vi] = _iog[vi];
}
fsas = new complex<double>[nsph]();
sscs = new double[nsph]();
sexs = new double[nsph]();
sabs = new double[nsph]();
sqscs = new double[nsph]();
sqexs = new double[nsph]();
sqabs = new double[nsph]();
gcsv = new double[nsph]();
rxx = new double[nsph]();
ryy = new double[nsph]();
rzz = new double[nsph]();
ros = new double[nsph]();
sas = new complex<double>**[nsph];
for (int si = 0; si < nsph; si++) {
sas[si] = new complex<double>*[2];
sas[si][0] = new complex<double>[2]();
sas[si][1] = new complex<double>[2]();
}
nlmmt = 2 * (l_max * (l_max + 2));
nsph = ns;
lm = l_max;
rmi = new complex<double>*[lm];
rei = new complex<double>*[lm];
for (int ri = 0; ri < lm; ri++) {
rmi[ri] = new complex<double>[nsph]();
rei[ri] = new complex<double>[nsph]();
}
w = new complex<double>*[nlmmt];
for (int wi = 0; wi < nlmmt; wi++) w[wi] = new complex<double>[4]();
int configurations = 0;
for (int ci = 1; ci <= nsph; ci++) {
if (_iog[ci - 1] >= ci) configurations++;
}
vints = new complex<double>*[nsph];
rc = new double*[configurations];
nshl = new int[configurations]();
iog = new int[nsph]();
int conf_index = 0;
for (int vi = 0; vi < nsph; vi++) {
vints[vi] = new complex<double>[16]();
iog[vi] = _iog[vi];
if (iog[vi] >= vi + 1) {
nshl[conf_index] = _nshl[conf_index];
rc[conf_index] = new double[_nshl[conf_index]]();
conf_index++;
}
}
fsas = new complex<double>[nsph]();
sscs = new double[nsph]();
sexs = new double[nsph]();
sabs = new double[nsph]();
sqscs = new double[nsph]();
sqexs = new double[nsph]();
sqabs = new double[nsph]();
gcsv = new double[nsph]();
rxx = new double[nsph]();
ryy = new double[nsph]();
rzz = new double[nsph]();
ros = new double[nsph]();
sas = new complex<double>**[nsph];
for (int si = 0; si < nsph; si++) {
sas[si] = new complex<double>*[2];
sas[si][0] = new complex<double>[2]();
sas[si][1] = new complex<double>[2]();
}
}
C1::~C1() {
delete[] rmi;
delete[] rei;
for (int wi = nlmmt - 1; wi > -1; wi--) delete[] w[wi];
for (int vi = nsph - 1; vi > - 1; vi--) {
delete[] rc[vi];
delete[] vints[vi];
}
delete[] rc;
delete[] vints;
for (int si = nsph - 1; si > -1; si--) {
delete[] sas[si][1];
delete[] sas[si][0];
delete[] sas[si];
}
delete[] sas;
delete[] fsas;
delete[] sscs;
delete[] sexs;
delete[] sabs;
delete[] sqscs;
delete[] sqexs;
delete[] sqabs;
delete[] gcsv;
delete[] rxx;
delete[] ryy;
delete[] rzz;
delete[] ros;
delete[] iog;
delete[] nshl;
delete[] rmi;
delete[] rei;
for (int wi = nlmmt - 1; wi > -1; wi--) delete[] w[wi];
int conf_index = 0;
for (int ci = 1; ci <= nsph; ci++) {
if (iog[ci] >= ci) {
delete[] rc[conf_index];
conf_index++;
}
}
delete[] rc;
for (int vi = nsph - 1; vi > - 1; vi--) {
delete[] vints[vi];
}
delete[] vints;
for (int si = nsph - 1; si > -1; si--) {
delete[] sas[si][1];
delete[] sas[si][0];
delete[] sas[si];
}
delete[] sas;
delete[] fsas;
delete[] sscs;
delete[] sexs;
delete[] sabs;
delete[] sqscs;
delete[] sqexs;
delete[] sqabs;
delete[] gcsv;
delete[] rxx;
delete[] ryy;
delete[] rzz;
delete[] ros;
delete[] iog;
delete[] nshl;
}
C1_AddOns::C1_AddOns(C4 *c4) {
nsph = c4->nsph;
lmpo = c4->lmpo;
nlemt = 2 * c4->nlem;
vh = new complex<double>[(nsph * nsph - 1) * c4->litpo]();
vj0 = new complex<double>[nsph * c4->lmtpo]();
vj = new complex<double>[1](); // QUESTION: is 1 really enough for a general case?
vyhj = new complex<double>[(nsph * nsph - 1) * c4->litpos]();
vyj0 = new complex<double>[nsph * c4->lmtpos]();
am0m = new complex<double>*[nlemt];
for (int ai = 0; ai < nlemt; ai++) {
am0m[ai] = new complex<double>[nlemt]();
}
vint = new complex<double>[16](); // QUESTION: is dimension 16 generally fixed?
vintm = new complex<double>[16]();
vintt = new complex<double>[16]();
vints = new complex<double>*[nsph];
for (int vi = 0; vi < nsph; vi++) vints[vi] = new complex<double>[16]();
fsac = new complex<double>*[2];
sac = new complex<double>*[2];
fsacm = new complex<double>*[2];
for (int fi = 0; fi < 2; fi++) {
fsac[fi] = new complex<double>[2]();
sac[fi] = new complex<double>[2]();
fsacm[fi] = new complex<double>[2]();
}
scscp = new complex<double>[2]();
ecscp = new complex<double>[2]();
scscpm = new complex<double>[2]();
ecscpm = new complex<double>[2]();
allocate_vectors(c4);
sscs = new double[nsph]();
ecscm = new double[2]();
scscm = new double[2]();
scsc = new double[2]();
ecsc = new double[2]();
nsph = c4->nsph;
lmpo = c4->lmpo;
nlemt = 2 * c4->nlem;
vh = new complex<double>[(nsph * nsph - 1) * c4->litpo]();
vj0 = new complex<double>[nsph * c4->lmtpo]();
vj = new complex<double>[1](); // QUESTION: is 1 really enough for a general case?
vyhj = new complex<double>[(nsph * nsph - 1) * c4->litpos]();
vyj0 = new complex<double>[nsph * c4->lmtpos]();
am0m = new complex<double>*[nlemt];
for (int ai = 0; ai < nlemt; ai++) {
am0m[ai] = new complex<double>[nlemt]();
}
vint = new complex<double>[16](); // QUESTION: is dimension 16 generally fixed?
vintm = new complex<double>[16]();
vintt = new complex<double>[16]();
vints = new complex<double>*[nsph];
for (int vi = 0; vi < nsph; vi++) vints[vi] = new complex<double>[16]();
fsac = new complex<double>*[2];
sac = new complex<double>*[2];
fsacm = new complex<double>*[2];
for (int fi = 0; fi < 2; fi++) {
fsac[fi] = new complex<double>[2]();
sac[fi] = new complex<double>[2]();
fsacm[fi] = new complex<double>[2]();
}
scscp = new complex<double>[2]();
ecscp = new complex<double>[2]();
scscpm = new complex<double>[2]();
ecscpm = new complex<double>[2]();
allocate_vectors(c4);
sscs = new double[nsph]();
ecscm = new double[2]();
scscm = new double[2]();
scsc = new double[2]();
ecsc = new double[2]();
}
C1_AddOns::~C1_AddOns() {
delete[] sscs;
delete[] vyj0;
delete[] vyhj;
for (int ii = lmpo - 1; ii > -1; ii--) delete[] ind3j[ii];
delete[] ind3j;
delete[] v3j0;
delete[] vh;
delete[] vj0;
for (int ai = nlemt - 1; ai > -1; ai--) {
delete[] am0m[ai];
}
delete am0m;
delete[] vint;
delete[] vintm;
delete[] vintt;
for (int vi = nsph - 1; vi > -1; vi--) delete[] vints[vi];
delete[] vints;
for (int fi = 1; fi > -1; fi--) {
delete[] fsac[fi];
delete[] sac[fi];
delete[] fsacm[fi];
}
delete[] fsac;
delete[] sac;
delete[] fsacm;
delete[] scscp;
delete[] ecscp;
delete[] scscpm;
delete[] ecscpm;
delete[] ecscm;
delete[] scscm;
delete[] scsc;
delete[] ecsc;
delete[] sscs;
delete[] vyj0;
delete[] vyhj;
for (int ii = lmpo - 1; ii > -1; ii--) delete[] ind3j[ii];
delete[] ind3j;
delete[] v3j0;
delete[] vh;
delete[] vj0;
for (int ai = nlemt - 1; ai > -1; ai--) {
delete[] am0m[ai];
}
delete am0m;
delete[] vint;
delete[] vintm;
delete[] vintt;
for (int vi = nsph - 1; vi > -1; vi--) delete[] vints[vi];
delete[] vints;
for (int fi = 1; fi > -1; fi--) {
delete[] fsac[fi];
delete[] sac[fi];
delete[] fsacm[fi];
}
delete[] fsac;
delete[] sac;
delete[] fsacm;
delete[] scscp;
delete[] ecscp;
delete[] scscpm;
delete[] ecscpm;
delete[] ecscm;
delete[] scscm;
delete[] scsc;
delete[] ecsc;
}
void C1_AddOns::allocate_vectors(C4 *c4) {
// This calculates the size of v3j0
int lm = (c4->li > c4->le) ? c4->li : c4->le;
const int nv3j = c4->nv3j;
v3j0 = new double[nv3j]();
ind3j = new int*[lmpo];
for (int ii = 0; ii < lmpo; ii++) ind3j[ii] = new int[c4->lm]();
// This calculates the size of vyhj
int ivy = (nsph * nsph - 1) * c4->litpos;
vyhj = new complex<double>[ivy]();
// This calculates the size of vyj0
ivy = c4->lmtpos * c4->nsph;
vyj0 = new complex<double>[ivy]();
// This calculates the size of v3j0
int lm = (c4->li > c4->le) ? c4->li : c4->le;
const int nv3j = c4->nv3j;
v3j0 = new double[nv3j]();
ind3j = new int*[lmpo];
for (int ii = 0; ii < lmpo; ii++) ind3j[ii] = new int[c4->lm]();
// This calculates the size of vyhj
int ivy = (nsph * nsph - 1) * c4->litpos;
vyhj = new complex<double>[ivy]();
// This calculates the size of vyj0
ivy = c4->lmtpos * c4->nsph;
vyj0 = new complex<double>[ivy]();
}
C2::C2(int ns, int nl, int npnt, int npntts) {
nsph = ns;
int max_n = (npnt > npntts) ? npnt : npntts;
nhspo = 2 * max_n - 1;
ris = new complex<double>[nhspo]();
dlri = new complex<double>[nhspo]();
vkt = new complex<double>[nsph]();
dc0 = new complex<double>[nl]();
vsz = new double[nsph]();
nsph = ns;
int max_n = (npnt > npntts) ? npnt : npntts;
nhspo = 2 * max_n - 1;
ris = new complex<double>[nhspo]();
dlri = new complex<double>[nhspo]();
vkt = new complex<double>[nsph]();
dc0 = new complex<double>[nl]();
vsz = new double[nsph]();
}
C2::~C2() {
delete[] ris;
delete[] dlri;
delete[] vkt;
delete[] dc0;
delete[] vsz;
delete[] ris;
delete[] dlri;
delete[] vkt;
delete[] dc0;
delete[] vsz;
}
C3::C3() {
tsas = new complex<double>*[2];
tsas[0] = new complex<double>[2];
tsas[1] = new complex<double>[2];
tfsas = complex<double>(0.0, 0.0);
gcs = 0.0;
scs = 0.0;
ecs = 0.0;
acs = 0.0;
tsas = new complex<double>*[2];
tsas[0] = new complex<double>[2];
tsas[1] = new complex<double>[2];
tfsas = complex<double>(0.0, 0.0);
gcs = 0.0;
scs = 0.0;
ecs = 0.0;
acs = 0.0;
}
C3::~C3() {
delete[] tsas[1];
delete[] tsas[0];
delete[] tsas;
delete[] tsas[1];
delete[] tsas[0];
delete[] tsas;
}
C6::C6(int lmtpo) {
rac3j = new double[lmtpo]();
rac3j = new double[lmtpo]();
}
C6::~C6() {
delete[] rac3j;
delete[] rac3j;
}
C9::C9(int ndi, int nlem, int ndit, int nlemt) {
gis = new complex<double>*[ndi];
gls = new complex<double>*[ndi];
for (int gi = 0; gi < ndi; gi++) {
gis[gi] = new complex<double>[nlem]();
gls[gi] = new complex<double>[nlem]();
}
sam = new complex<double>*[ndit];
for (int si = 0; si < ndit; si++) sam[si] = new complex<double>[nlemt]();
gis_size_0 = ndi;
sam_size_0 = ndit;
gis = new complex<double>*[ndi];
gls = new complex<double>*[ndi];
for (int gi = 0; gi < ndi; gi++) {
gis[gi] = new complex<double>[nlem]();
gls[gi] = new complex<double>[nlem]();
}
sam = new complex<double>*[ndit];
for (int si = 0; si < ndit; si++) sam[si] = new complex<double>[nlemt]();
gis_size_0 = ndi;
sam_size_0 = ndit;
}
C9::~C9() {
for (int gi = gis_size_0 - 1; gi > -1; gi--) {
delete[] gis[gi];
delete[] gls[gi];
}
delete[] gis;
delete[] gls;
for (int si = sam_size_0 - 1; si > -1; si--) delete[] sam[si];
delete[] sam;
for (int gi = gis_size_0 - 1; gi > -1; gi--) {
delete[] gis[gi];
delete[] gls[gi];
}
delete[] gis;
delete[] gls;
for (int si = sam_size_0 - 1; si > -1; si--) delete[] sam[si];
delete[] sam;
}
src/libnptm/Configuration.cpp
+ 6
− 0
View file @ 1aaea3d6
......@@ -1067,6 +1067,12 @@ void ScattererConfiguration::write_formatted(string file_name) {
);
break;
}
// Clean memory
delete[] xi_vec;
delete[] pu_vec;
delete[] ev_vec;
delete[] wn_vec;
delete[] wl_vec;
} else { // idfc < 0, Dielectric functions are at XIP and XI is scale for dimensions
double pu, wn;
xi_vec = scale_vec;
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment