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

Transform the AM matrix into a contiguous memory object and skip copying before LAPCK call

parent 8295d147
No related branches found
No related tags found
No related merge requests found
...@@ -144,9 +144,11 @@ void cluster(string config_file, string data_file, string output_path) { ...@@ -144,9 +144,11 @@ void cluster(string config_file, string data_file, string output_path) {
if (sconf->iog_vec[ci -1] >= ci) configurations++; if (sconf->iog_vec[ci -1] >= ci) configurations++;
} }
C2 *c2 = new C2(nsph, configurations, npnt, npntts); C2 *c2 = new C2(nsph, configurations, npnt, npntts);
dcomplex **am = new dcomplex*[mxndm]; lapack_int ndit = 2 * nsph * c4->nlim;
for (int ai = 0; ai < mxndm; ai++) { dcomplex *am_vector = new dcomplex[ndit * ndit]();
am[ai] = new dcomplex[mxndm](); dcomplex **am = new dcomplex*[ndit];
for (int ai = 0; ai < ndit; ai++) {
am[ai] = (am_vector + ai * ndit);
} }
const int ndi = c4->nsph * c4->nlim; const int ndi = c4->nsph * c4->nlim;
C9 *c9 = new C9(ndi, c4->nlem, 2 * ndi, 2 * c4->nlem); C9 *c9 = new C9(ndi, c4->nlem, 2 * ndi, 2 * c4->nlem);
...@@ -359,7 +361,6 @@ void cluster(string config_file, string data_file, string output_path) { ...@@ -359,7 +361,6 @@ void cluster(string config_file, string data_file, string output_path) {
if (jer != 0) break; if (jer != 0) break;
} // i132 loop } // i132 loop
cms(am, c1, c1ao, c4, c6); cms(am, c1, c1ao, c4, c6);
lapack_int ndit = 2 * nsph * c4->nlim;
invert_matrix(am, ndit, jer, mxndm); invert_matrix(am, ndit, jer, mxndm);
if (jer != 0) break; // jxi488 loop: goes to memory clean if (jer != 0) break; // jxi488 loop: goes to memory clean
ztm(am, c1, c1ao, c4, c6, c9); ztm(am, c1, c1ao, c4, c6, c9);
...@@ -937,10 +938,7 @@ void cluster(string config_file, string data_file, string output_path) { ...@@ -937,10 +938,7 @@ void cluster(string config_file, string data_file, string output_path) {
delete[] zpv[zi]; delete[] zpv[zi];
} }
delete[] zpv; delete[] zpv;
for (int ai = mxndm - 1; ai > -1; ai--) { delete[] am_vector;
delete[] am[ai];
//delete[] tam[ai];
}
delete[] am; delete[] am;
//delete[] tam; //delete[] tam;
delete[] gaps; delete[] gaps;
......
...@@ -23,18 +23,12 @@ ...@@ -23,18 +23,12 @@
#ifdef USE_LAPACK #ifdef USE_LAPACK
void zinvert(dcomplex **mat, lapack_int n, int &jer) { void zinvert(dcomplex **mat, lapack_int n, int &jer) {
jer = 0; jer = 0;
dcomplex *arr = new dcomplex[n * n]; dcomplex *arr = &(mat[0][0]);
const dcomplex uim = 0.0 + 1.0 * I; const dcomplex uim = 0.0 + 1.0 * I;
#ifdef USE_MKL #ifdef USE_MKL
MKL_Complex16 *arr2 = (MKL_Complex16 *) arr; MKL_Complex16 *arr2 = (MKL_Complex16 *) arr;
#endif #endif
for (lapack_int i = 0; i < n; i++) {
for (lapack_int j = 0; j < n; j++) {
lapack_int idx = i + n * j;
arr[idx] = mat[j][i];
}
}
lapack_int* IPIV = new lapack_int[n](); lapack_int* IPIV = new lapack_int[n]();
...@@ -46,13 +40,6 @@ void zinvert(dcomplex **mat, lapack_int n, int &jer) { ...@@ -46,13 +40,6 @@ void zinvert(dcomplex **mat, lapack_int n, int &jer) {
LAPACKE_zgetri(LAPACK_ROW_MAJOR, n, arr, n, IPIV); LAPACKE_zgetri(LAPACK_ROW_MAJOR, n, arr, n, IPIV);
#endif #endif
for (lapack_int i = 0; i < n; i++) {
for (lapack_int j = 0; j < n; j++) {
lapack_int idx = i + n * j;
mat[j][i] = arr[idx];
}
}
delete[] IPIV; delete[] IPIV;
delete[] arr;
} }
#endif #endif
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment