Skip to content
GitLab
Explore
Sign in
Primary navigation
Search or go to…
Project
E
euroexa_GaiaGsrParSim
Manage
Activity
Members
Labels
Plan
Issues
Issue boards
Milestones
Wiki
Code
Merge requests
Repository
Branches
Commits
Tags
Repository graph
Compare revisions
Snippets
Build
Pipelines
Jobs
Pipeline schedules
Artifacts
Deploy
Releases
Container registry
Model registry
Operate
Environments
Monitor
Incidents
Analyze
Value stream analytics
Contributor analytics
CI/CD analytics
Repository analytics
Model experiments
Help
Help
Support
GitLab documentation
Compare GitLab plans
Community forum
Contribute to GitLab
Provide feedback
Keyboard shortcuts
?
Snippets
Groups
Projects
Show more breadcrumbs
Fabio Roberto Vitello
euroexa_GaiaGsrParSim
Commits
378fbe40
Commit
378fbe40
authored
6 years ago
by
Fabio Roberto Vitello
Browse files
Options
Downloads
Patches
Plain Diff
Divisione in task
parent
692132e7
No related branches found
No related tags found
No related merge requests found
Changes
3
Expand all
Show whitespace changes
Inline
Side-by-side
Showing
3 changed files
aprod.c
+693
-422
693 additions, 422 deletions
aprod.c
lsqr.c
+1314
-1197
1314 additions, 1197 deletions
lsqr.c
solvergaiaSim.c
+2575
-2194
2575 additions, 2194 deletions
solvergaiaSim.c
with
4582 additions
and
3813 deletions
aprod.c
+
693
−
422
View file @
378fbe40
...
...
@@ -12,7 +12,7 @@ void aprod(int mode, long int m, long int n, double *vVect, double *knownTerms,
// Parallel definitions
int
myid
,
nproc
;
long
int
*
mapNoss
,
*
mapNcoeff
;
int
nthreads
,
tid
;
int
nthreads
,
tid
,
ntasks
;
long
**
mapForThread
;
/// struct comData *comlsqr;
//
...
...
@@ -38,6 +38,8 @@ nproc=comlsqr.nproc;
mapNcoeff
=
comlsqr
.
mapNcoeff
;
mapNoss
=
comlsqr
.
mapNoss
;
nthreads
=
comlsqr
.
nthreads
;
ntasks
=
comlsqr
.
ntasks
;
mapForThread
=
comlsqr
.
mapForThread
;
int
multMI
=
comlsqr
.
multMI
;
long
nparam
=
comlsqr
.
parOss
;
...
...
@@ -68,14 +70,12 @@ setBound[0]=comlsqr.setBound[0];
long
offsetInstrParam
=
comlsqr
.
offsetInstrParam
;
long
offsetGlobParam
=
comlsqr
.
offsetGlobParam
;
nthreads
=
1
;
//
nthreads
=
1;
tid
=
0
;
FILE
*
fp1
,
*
fp2
;
// fp1=fopen("test1_aprod","w");
// fp2=fopen("test2_aprod","w");
if
(
mode
!=
1
&&
mode
!=
2
)
{
printf
(
"ERROR: Invalid mode=%d in aprod function
\n
"
,
mode
);
...
...
@@ -87,16 +87,18 @@ myid=comlsqr.myid;
if
(
mode
==
1
)
{
time_t
startTime
=
time
(
NULL
);
#pragma omp parallel private(myid, sum,
k, l1, l2, l, j,tid,nthreads,i2,na) shared(mapNoss,instrCol,comlsqr,vVect,systemMatrix,matrixIndex,knownTerms,j2)
////
#pragma omp parallel private(myid, sum, k, l1, l2, l, j,
tid,
nthreads,
i2,
na) shared(mapNoss,
instrCol,
comlsqr,
vVect,
systemMatrix,
matrixIndex,
knownTerms,
j2)
{
myid
=
comlsqr
.
myid
;
/*
//FV_ EDIT ompSs
if (comlsqr.itn == 1)
{
#ifdef OMP
tid = omp_get_thread_num();
nthreads = omp_get_num_threads();
#endif
}
}*/
if
(
comlsqr
.
itn
==
1
&&
debugMode
)
printf
(
"PE=%d Aprod1 OpenMP num of threads =%d from thread =%d icycle=%ld comlsqr.itn=%d
\n
"
,
myid
,
nthreads
,
tid
,
i
,
comlsqr
.
itn
);
long
miValAstro
=
0
;
...
...
@@ -117,31 +119,41 @@ if(mode==1)
int
nGlobVal
=
nAstroPSolved
+
nAttP
+
nInstrPSolved
;
jstartAstro
=
miValAstro
-
offLocalAstro
;
//FV_ EDIT ompSs
#pragma omp for
for
(
long
ix
=
0
;
ix
<
mapNoss
[
myid
];
ix
++
){
for
(
int
nt
=
0
;
nt
<
ntasks
;
nt
++
)
{
#pragma omp task //shared(matrixIndex,mapForThread,vVect,systemMatrix,knownTerms) in(nt,multMI,offLocalAstro, lset)
{
for
(
long
ix
=
mapForThread
[
nt
][
0
];
ix
<
mapForThread
[
nt
][
2
];
ix
++
)
{
sum
=
0
.;
/////////////////////////////////////////////////////
/// Mode 1 Astrometric Sect
if
(
nAstroPSolved
){
if
(
nAstroPSolved
)
{
lset
=
ix
*
nparam
;
if
(
matrixIndex
[
multMI
*
ix
]
!=
miValAstro
){
if
(
matrixIndex
[
multMI
*
ix
]
!=
miValAstro
)
{
miValAstro
=
matrixIndex
[
multMI
*
ix
];
jstartAstro
=
miValAstro
-
offLocalAstro
;
}
for
(
long
jx
=
jstartAstro
;
jx
<
jstartAstro
+
nAstroPSolved
;
jx
++
){
for
(
long
jx
=
jstartAstro
;
jx
<
jstartAstro
+
nAstroPSolved
;
jx
++
)
{
sum
=
sum
+
systemMatrix
[
lset
]
*
vVect
[
jx
];
lset
++
;
}
}
//////////////////////////////////////////////////////
/// Mode 1 Attitude Sect
if
(
nAttP
)
{
if
(
nAttP
)
{
lset
=
ix
*
nparam
+
nAstroPSolved
;
miValAtt
=
matrixIndex
[
multMI
*
ix
+
(
multMI
-
1
)];
for
(
int
nax
=
0
;
nax
<
nAttAxes
;
nax
++
){
for
(
int
nax
=
0
;
nax
<
nAttAxes
;
nax
++
)
{
jstartAtt
=
miValAtt
+
offLocalAtt
+
nax
*
nDegFreedomAtt
;
for
(
long
inpax
=
jstartAtt
;
inpax
<
jstartAtt
+
nAttParAxis
;
inpax
++
)
{
...
...
@@ -152,11 +164,13 @@ if (nAttP){
}
//////////////////////////////////////////////////////
/// Mode 1 Instrument Sect
if
(
nInstrPSolved
){
if
(
nInstrPSolved
)
{
lset
=
ix
*
nparam
+
nInstrVal
;
long
iiVal
=
ix
*
nInstrPSolved
;
for
(
int
inInstr
=
0
;
inInstr
<
nInstrPSolved
;
inInstr
++
){
for
(
int
inInstr
=
0
;
inInstr
<
nInstrPSolved
;
inInstr
++
)
{
ixInstr
=
offLocalInstr
+
instrCol
[
iiVal
+
inInstr
];
sum
=
sum
+
systemMatrix
[
lset
]
*
vVect
[
ixInstr
];
lset
++
;
...
...
@@ -164,61 +178,155 @@ if (nInstrPSolved){
}
//////////////////////////////////////////////////////
/// Mode 1 Global sect
if
(
nGlobP
){
if
(
nGlobP
)
{
lset
=
ix
*
nparam
+
nGlobVal
;
for
(
long
inGlob
=
offLocalGlob
;
inGlob
<
offLocalGlob
+
nGlobP
;
inGlob
++
){
for
(
long
inGlob
=
offLocalGlob
;
inGlob
<
offLocalGlob
+
nGlobP
;
inGlob
++
)
{
sum
=
sum
+
systemMatrix
[
lset
]
*
vVect
[
inGlob
];
lset
++
;
}
}
//////////////////////////////////////////////////////
knownTerms
[
ix
]
+=
sum
;
}
//for ix
}
#pragma omp taskwait
}
/*
//#pragma omp for
for (long ix = 0; ix < mapNoss[myid]; ix++)
{
sum = 0.;
/////////////////////////////////////////////////////
/// Mode 1 Astrometric Sect
if (nAstroPSolved)
{
lset = ix * nparam;
if (matrixIndex[multMI * ix] != miValAstro)
{
miValAstro = matrixIndex[multMI * ix];
jstartAstro = miValAstro - offLocalAstro;
}
for (long jx = jstartAstro; jx < jstartAstro + nAstroPSolved; jx++)
{
sum = sum + systemMatrix[lset] * vVect[jx];
lset++;
}
}
//////////////////////////////////////////////////////
/// Mode 1 Attitude Sect
if (nAttP)
{
lset = ix * nparam + nAstroPSolved;
miValAtt = matrixIndex[multMI * ix + (multMI - 1)];
for (int nax = 0; nax < nAttAxes; nax++)
{
jstartAtt = miValAtt + offLocalAtt + nax * nDegFreedomAtt;
for (long inpax = jstartAtt; inpax < jstartAtt + nAttParAxis; inpax++)
{
sum = sum + systemMatrix[lset] * vVect[inpax];
lset++;
}
}
}
//////////////////////////////////////////////////////
/// Mode 1 Instrument Sect
if (nInstrPSolved)
{
lset = ix * nparam + nInstrVal;
long iiVal = ix * nInstrPSolved;
for (int inInstr = 0; inInstr < nInstrPSolved; inInstr++)
{
ixInstr = offLocalInstr + instrCol[iiVal + inInstr];
sum = sum + systemMatrix[lset] * vVect[ixInstr];
lset++;
}
}
//////////////////////////////////////////////////////
/// Mode 1 Global sect
if (nGlobP)
{
lset = ix * nparam + nGlobVal;
for (long inGlob = offLocalGlob; inGlob < offLocalGlob + nGlobP; inGlob++)
{
sum = sum + systemMatrix[lset] * vVect[inGlob];
lset++;
}
}
//////////////////////////////////////////////////////
knownTerms[ix] += sum;
} //for ix
*/
/// Mode 1 ExtConstr
if
(
nEqExtConstr
){
if
(
nEqExtConstr
)
{
long
offExtAtt
;
long
offExtAttConstr
=
VrIdAstroPDimMax
*
nAstroPSolved
+
startingAttColExtConstr
;
long
vVIx
;
long
ktIx
=
mapNoss
[
myid
];
long
offExtConstr
;
for
(
int
iexc
=
0
;
iexc
<
nEqExtConstr
;
iexc
++
){
for
(
int
iexc
=
0
;
iexc
<
nEqExtConstr
;
iexc
++
)
{
sum
=
0
.
0
;
offExtConstr
=
mapNcoeff
[
myid
]
+
iexc
*
nOfElextObs
;
#pragma omp for
//FV_ EDIT ompSs
//#pragma omp for
for
(
int
j3
=
0
;
j3
<
numOfExtStar
*
nAstroPSolved
;
j3
++
)
sum
+=
systemMatrix
[
offExtConstr
+
j3
]
*
vVect
[
j3
];
for
(
int
nax
=
0
;
nax
<
nAttAxes
;
nax
++
){
for
(
int
nax
=
0
;
nax
<
nAttAxes
;
nax
++
)
{
offExtAtt
=
offExtConstr
+
numOfExtStar
*
nAstroPSolved
+
nax
*
numOfExtAttCol
;
vVIx
=
offExtAttConstr
+
nax
*
nDegFreedomAtt
;
#pragma omp for
//FV_ EDIT ompSs
//#pragma omp for
for
(
int
j3
=
0
;
j3
<
numOfExtAttCol
;
j3
++
)
sum
+=
systemMatrix
[
offExtAtt
+
j3
]
*
vVect
[
vVIx
+
j3
];
}
#pragma omp atomic
//FV_ EDIT ompSs
//#pragma omp atomic
knownTerms
[
ktIx
+
iexc
]
+=
sum
;
}
//for iexc
}
//////////////////////////////////////////////////////
/// Mode 1 BarConstr
if
(
nEqBarConstr
){
if
(
nEqBarConstr
)
{
long
offBarConstr
=
mapNcoeff
[
myid
]
+
nOfElextObs
*
nEqExtConstr
;
long
offBarConstrIx
;
long
ktIx
=
mapNoss
[
myid
]
+
nEqExtConstr
;
for
(
int
iexc
=
0
;
iexc
<
nEqBarConstr
;
iexc
++
){
for
(
int
iexc
=
0
;
iexc
<
nEqBarConstr
;
iexc
++
)
{
sum
=
0
.
0
;
offBarConstrIx
=
offBarConstr
+
iexc
*
nOfElBarObs
;
#pragma omp for
//FV_ EDIT ompSs
//#pragma omp for
for
(
int
j3
=
0
;
j3
<
numOfBarStar
*
nAstroPSolved
;
j3
++
)
sum
+=
systemMatrix
[
offBarConstrIx
+
j3
]
*
vVect
[
j3
];
#pragma omp atomic
//FV_ EDIT ompSs
//#pragma omp atomic
knownTerms
[
ktIx
+
iexc
]
+=
sum
;
}
//for iexc
}
//////////////////////////////////////////////////////
/// Mode 1 InstrConstr
if
(
nOfInstrConstr
){
if
(
nOfInstrConstr
)
{
long
offSetInstr
=
0
;
long
offSetInstrInc
=
0
;
long
offSetInstrInc1
=
0
;
...
...
@@ -227,77 +335,164 @@ if(nOfInstrConstr){
long
offSetInstrConstr
=
mapNcoeff
[
myid
]
+
nOfElextObs
*
nEqExtConstr
+
nOfElBarObs
*
nEqBarConstr
;
long
offSetInstrConstr1
=
VrIdAstroPDimMax
*
nAstroPSolved
+
nDegFreedomAtt
*
nAttAxes
;
long
ktIx
=
mapNoss
[
myid
]
+
nEqExtConstr
+
nEqBarConstr
;
for
(
int
i1
=
myid
;
i1
<
nOfInstrConstr
;
i1
=
i1
+
nproc
){
for
(
int
i1
=
myid
;
i1
<
nOfInstrConstr
;
i1
=
i1
+
nproc
)
{
sum
=
0
.
0
;
offSetInstrInc
=
offSetInstrConstr
;
offSetInstr
=
0
;
for
(
int
m
=
0
;
m
<
i1
;
m
++
){
for
(
int
m
=
0
;
m
<
i1
;
m
++
)
{
offSetInstrInc
+=
instrConstrIlung
[
m
];
offSetInstr
+=
instrConstrIlung
[
m
];
}
offvV
=
mapNoss
[
myid
]
*
nInstrPSolved
+
offSetInstr
;
#pragma omp for
for
(
int
j3
=
0
;
j3
<
instrConstrIlung
[
i1
];
j3
++
){
//FV_ EDIT ompSs
//#pragma omp for
for
(
int
j3
=
0
;
j3
<
instrConstrIlung
[
i1
];
j3
++
)
{
vVix
=
instrCol
[
offvV
+
j3
];
sum
+=
systemMatrix
[
offSetInstrInc
+
j3
]
*
vVect
[
offSetInstrConstr1
+
vVix
];
}
#pragma omp atomic
//FV_ EDIT ompSs
//#pragma omp atomic
knownTerms
[
ktIx
+
i1
]
+=
sum
;
}
}
//////////////////////////////////////////////////////
}
//pragma
*
ompSec
+=
time
(
NULL
)
-
startTime
;
if
(
comlsqr
.
itn
<=
2
&&
(
myid
==
0
||
myid
==
nproc
-
1
||
debugMode
==
1
))
printf
(
"PE=%d AprodTiming: mode=1 OmpSec=%ld
\n
"
,
myid
,
time
(
NULL
)
-
startTime
);
if
(
comlsqr
.
itn
<=
2
&&
(
myid
==
0
||
myid
==
nproc
-
1
||
debugMode
==
1
))
printf
(
"PE=%d AprodTiming: mode=1 OmpSec=%ld
\n
"
,
myid
,
time
(
NULL
)
-
startTime
);
}
else
//mode==2
{
//if(mode
time_t
startTime
=
time
(
NULL
);
#pragma omp parallel private(myid,yi, sum, k, l1, l2, l, j,localSum,i,tid,nthreads) shared(mapNoss,instrCol,comlsqr,vVect,systemMatrix,matrixIndex,knownTerms)
//FV_ EDIT ompSs
// #pragma omp parallel private(myid, yi, sum, k, l1, l2, l, j, localSum, i, tid, nthreads) shared(mapNoss, instrCol, comlsqr, vVect, systemMatrix, matrixIndex, knownTerms)
{
int
count
=
0
;
myid
=
comlsqr
.
myid
;
//qui lavorare
/*
#ifdef OMP
tid = omp_get_thread_num();
nthreads = omp_get_num_threads();
#endif
*/
/////////////////////////////////////////////////////
/// Mode 2 Astrometric Sect
if
(
nAstroPSolved
){
if
(
nAstroPSolved
)
{
long
offLocalAstro
=
comlsqr
.
mapStar
[
myid
][
0
]
*
nAstroPSolved
;
// for (long ix = mapForThread[tid][0]; ix < mapForThread[tid][2]; ix++)
long
lset
;
for
(
int
nt
=
0
;
nt
<
ntasks
;
nt
++
)
{
lset
=
mapForThread
[
nt
][
0
]
*
nparam
;
#pragma omp task shared(matrixIndex,mapForThread,vVect) in(nt,multMI,offLocalAstro, lset)
{
double
taskLocalSum
;
long
jstartAstro
=
0
;
long
lset
=
mapForThread
[
tid
][
0
]
*
nparam
;
for
(
long
ix
=
mapForThread
[
tid
][
0
];
ix
<
mapForThread
[
tid
][
2
];
ix
++
){
for
(
long
ix
=
mapForThread
[
nt
][
0
];
ix
<
mapForThread
[
nt
][
2
];
ix
++
)
{
long
miValAstro
=
0
;
if
(
matrixIndex
[
multMI
*
ix
]
!=
miValAstro
){
if
(
matrixIndex
[
multMI
*
ix
]
!=
miValAstro
)
{
miValAstro
=
matrixIndex
[
multMI
*
ix
];
jstartAstro
=
miValAstro
-
offLocalAstro
;
}
for
(
long
jx
=
jstartAstro
;
jx
<
jstartAstro
+
nAstroPSolved
;
jx
++
){
for
(
long
jx
=
jstartAstro
;
jx
<
jstartAstro
+
nAstroPSolved
;
jx
++
)
{
taskLocalSum
=
systemMatrix
[
lset
]
*
knownTerms
[
ix
];
vVect
[
jx
]
+=
taskLocalSum
;
lset
++
;
}
//for jx
lset
+=
nparam
-
nAstroPSolved
;
}
//for ix
}
#pragma omp taskwait
}
/*
for (long ix = 0; ix < comlsqr.mapNoss[myid]; ix++)
{
long miValAstro = 0;
if (matrixIndex[multMI * ix] != miValAstro)
{
miValAstro = matrixIndex[multMI * ix];
jstartAstro = miValAstro - offLocalAstro;
}
for (long jx = jstartAstro; jx < jstartAstro + nAstroPSolved; jx++)
{
localSum = systemMatrix[lset] * knownTerms[ix];
vVect[jx] += localSum;
lset++;
} //for jx
lset += nparam - nAstroPSolved;
} //for ix
*/
}
/////////////////////////////////////////////////////
}
//pragma
/////////////////////////////////////////////////////
/// Mode 2 Attitude Sect
if
(
nAttP
){
long
lset
;
if
(
nAttP
)
{
long
mix
;
long
offj
=
(
localAstroMax
-
offsetAttParam
);
long
vVix
;
for
(
int
nt
=
0
;
nt
<
ntasks
;
nt
++
)
{
long
lset
;
// #pragma omp task shared(nparam,nInstrPSolved,matrixIndex) in (nt) commutative(vVect)
#pragma omp task commutative(vVect)
{
for
(
long
ix
=
mapForThread
[
nt
][
0
];
ix
<
mapForThread
[
nt
][
2
];
ix
++
)
{
lset
=
nparam
*
ix
+
nAstroPSolved
;
mix
=
matrixIndex
[
multMI
*
ix
+
(
multMI
-
1
)]
+
offj
;
for
(
int
ly
=
0
;
ly
<
nAttAxes
;
ly
++
)
{
vVix
=
mix
+
ly
*
nDegFreedomAtt
;
for
(
int
lx
=
0
;
lx
<
nAttParAxis
;
lx
++
)
{
localSum
=
systemMatrix
[
lset
]
*
knownTerms
[
ix
];
///#pragma omp atomic
vVect
[
vVix
]
+=
localSum
;
lset
++
;
vVix
++
;
}
//for lx
}
//for ly
}
//for ix
}
#pragma omp taskwait
}
/*
///#pragma omp for
for
(
long
ix
=
0
;
ix
<
mapNoss
[
myid
];
ix
++
){
for (long ix = 0; ix < mapNoss[myid]; ix++)
{
lset = nparam * ix + nAstroPSolved;
mix = matrixIndex[multMI * ix + (multMI - 1)] + offj;
for
(
int
ly
=
0
;
ly
<
nAttAxes
;
ly
++
){
for (int ly = 0; ly < nAttAxes; ly++)
{
vVix = mix + ly * nDegFreedomAtt;
for
(
int
lx
=
0
;
lx
<
nAttParAxis
;
lx
++
){
for (int lx = 0; lx < nAttParAxis; lx++)
{
localSum = systemMatrix[lset] * knownTerms[ix];
///#pragma omp atomic
vVect[vVix] += localSum;
...
...
@@ -306,20 +501,49 @@ if (nAttP){
} //for lx
} //for ly
} //for ix
*/
}
/////////////////////////////////////////////////////
/// Mode 2 Instrument Sect
if
(
nOfInstrConstr
){
long
lset
;
if
(
nOfInstrConstr
)
{
int
offLset
=
nAstroPSolved
+
nAttP
;
long
iix
;
long
offj
=
offsetInstrParam
+
(
localAstroMax
-
offsetAttParam
);
long
vVix
;
for
(
int
nt
=
0
;
nt
<
ntasks
;
nt
++
)
{
// #pragma omp task shared(vVect,nparam,offLset,nInstrPSolved) in (nt)
#pragma omp task commutative(vVect)
{
long
lset
;
for
(
long
ix
=
mapForThread
[
nt
][
0
];
ix
<
mapForThread
[
nt
][
2
];
ix
++
)
{
lset
=
nparam
*
ix
+
offLset
;
iix
=
ix
*
nInstrPSolved
;
for
(
int
ly
=
0
;
ly
<
nInstrPSolved
;
ly
++
)
{
vVix
=
offj
+
instrCol
[
iix
+
ly
];
localSum
=
systemMatrix
[
lset
]
*
knownTerms
[
ix
];
///#pragma omp atomic
vVect
[
vVix
]
+=
localSum
;
lset
++
;
}
//for ly
}
//for ix
}
#pragma omp taskwait
}
/*
///#pragma omp for
for
(
long
ix
=
0
;
ix
<
mapNoss
[
myid
];
ix
++
){
for (long ix = 0; ix < mapNoss[myid]; ix++)
{
lset = nparam * ix + offLset;
iix = ix * nInstrPSolved;
for
(
int
ly
=
0
;
ly
<
nInstrPSolved
;
ly
++
){
for (int ly = 0; ly < nInstrPSolved; ly++)
{
vVix = offj + instrCol[iix + ly];
localSum = systemMatrix[lset] * knownTerms[ix];
///#pragma omp atomic
...
...
@@ -327,103 +551,144 @@ if(nOfInstrConstr){
lset++;
} //for ly
} //for ix
*/
}
/////////////////////////////////////////////////////
///} //pragma
*
ompSec
+=
time
(
NULL
)
-
startTime
;
if
(
comlsqr
.
itn
<=
2
&&
(
myid
==
0
||
myid
==
nproc
-
1
||
debugMode
==
1
)
)
printf
(
"PE=%d AprodTiming: mode=2.1 OmpSec=%ld
\n
"
,
myid
,
time
(
NULL
)
-
startTime
);
if
(
comlsqr
.
itn
<=
2
&&
(
myid
==
0
||
myid
==
nproc
-
1
||
debugMode
==
1
))
printf
(
"PE=%d AprodTiming: mode=2.1 OmpSec=%ld
\n
"
,
myid
,
time
(
NULL
)
-
startTime
);
/////////////////////////////////////////////////////
/// Mode 2 Global Sect
if
(
nGlobP
){
if
(
nGlobP
)
{
time_t
startTime_nGlobP
=
time
(
NULL
);
long
lset
;
int
offLset
=
nAstroPSolved
+
nAttP
+
nInstrPSolved
;
long
offj
=
offsetGlobParam
+
(
localAstroMax
-
offsetAttParam
);
long
vVix
;
for
(
long
ix
=
0
;
ix
<
mapNoss
[
myid
];
ix
++
){
for
(
long
ix
=
0
;
ix
<
mapNoss
[
myid
];
ix
++
)
{
lset
=
nparam
*
ix
+
offLset
;
for
(
long
ly
=
0
;
ly
<
nGlobP
;
ly
++
){
for
(
long
ly
=
0
;
ly
<
nGlobP
;
ly
++
)
{
vVix
=
offj
+
ly
;
localSum
=
systemMatrix
[
lset
]
*
knownTerms
[
ix
];
vVect
[
vVix
]
+=
localSum
;
lset
++
;
}
//for ly
}
//for ix
if
(
comlsqr
.
itn
<=
2
&&
(
myid
==
0
||
myid
==
nproc
-
1
||
debugMode
==
1
))
printf
(
"
\t
PE=%d AprodTiming: nGlobP sec=%ld
\n
"
,
myid
,
time
(
NULL
)
-
startTime_nGlobP
);
}
////////////////////////////////////////////////////
startTime
=
time
(
NULL
);
#pragma omp parallel private(myid,yi,localSum,tid,nthreads,i2,j2,na) shared(mapNoss,vVect,systemMatrix,knownTerms,k,j3)
//FV_ EDIT ompSs
//
#pragma omp parallel private(myid,
yi,
localSum,
tid,
nthreads,
i2,
j2,
na) shared(mapNoss,
vVect,
systemMatrix,
knownTerms,
k,
j3)
{
myid
=
comlsqr
.
myid
;
#ifdef OMP
//FV_ EDIT ompSs
/*
//#ifdef OMP
tid = omp_get_thread_num();
nthreads = omp_get_num_threads();
#endif
*/
localSum
=
0
.
0
;
//////////////////////////////////////////////////////
/// Mode 2 ExtConstr
if
(
nEqExtConstr
){
if
(
nEqExtConstr
)
{
time_t
startTime_nEqExtConstr
=
time
(
NULL
);
long
offExtStarConstrEq
;
long
off2
;
long
off3
;
for
(
int
ix
=
0
;
ix
<
nEqExtConstr
;
ix
++
){
//stars
for
(
int
ix
=
0
;
ix
<
nEqExtConstr
;
ix
++
)
{
//stars
yi
=
knownTerms
[
mapNoss
[
myid
]
+
ix
];
offExtStarConstrEq
=
mapNcoeff
[
myid
]
+
ix
*
nOfElextObs
;
//Star offset on the row ix (all the row)
#pragma omp for
for
(
long
yx
=
0
;
yx
<
numOfExtStar
;
yx
++
){
//FV_ EDIT ompSs
//#pragma omp for
for
(
long
yx
=
0
;
yx
<
numOfExtStar
;
yx
++
)
{
off3
=
yx
*
nAstroPSolved
;
off2
=
offExtStarConstrEq
+
off3
;
for
(
int
j2
=
0
;
j2
<
nAstroPSolved
;
j2
++
){
for
(
int
j2
=
0
;
j2
<
nAstroPSolved
;
j2
++
)
{
localSum
=
systemMatrix
[
off2
+
j2
]
*
yi
;
vVect
[
j2
+
off3
]
+=
localSum
;
}
}
#pragma omp barrier
//FV_ EDIT ompSs
//#pragma omp barrier
}
//for ix
long
offExtAttConstrEq
;
long
offExtUnk
;
long
off1
=
VrIdAstroPDimMax
*
nAstroPSolved
+
startingAttColExtConstr
;
for
(
int
ix
=
0
;
ix
<
nEqExtConstr
;
ix
++
){
//att
for
(
int
ix
=
0
;
ix
<
nEqExtConstr
;
ix
++
)
{
//att
yi
=
knownTerms
[
mapNoss
[
myid
]
+
ix
];
offExtAttConstrEq
=
mapNcoeff
[
myid
]
+
ix
*
nOfElextObs
;
//Att offset on the row ix (all the row)
offExtAttConstrEq
+=
numOfExtStar
*
nAstroPSolved
;
//Att offset inside ix row
for
(
int
nax
=
0
;
nax
<
nAttAxes
;
nax
++
){
for
(
int
nax
=
0
;
nax
<
nAttAxes
;
nax
++
)
{
offExtUnk
=
off1
+
nax
*
nDegFreedomAtt
;
// Att offset for Unk array on extConstr
off2
=
offExtAttConstrEq
+
nax
*
numOfExtAttCol
;
#pragma omp for
for
(
int
jx
=
0
;
jx
<
numOfExtAttCol
;
jx
++
){
//FV_ EDIT ompSs
//#pragma omp for
for
(
int
jx
=
0
;
jx
<
numOfExtAttCol
;
jx
++
)
{
localSum
=
systemMatrix
[
off2
+
jx
]
*
yi
;
vVect
[
offExtUnk
+
jx
]
+=
localSum
;
}
#pragma omp barrier
//FV_ EDIT ompSs
//#pragma omp barrier
}
}
if
(
comlsqr
.
itn
<=
2
&&
(
myid
==
0
||
myid
==
nproc
-
1
||
debugMode
==
1
))
printf
(
"
\t
PE=%d AprodTiming: nEqExtConstr sec=%ld
\n
"
,
myid
,
time
(
NULL
)
-
startTime_nEqExtConstr
);
}
//////////////////////////////////////////////////////
/// Mode 2 BarConstr
if
(
nEqBarConstr
){
if
(
nEqBarConstr
)
{
time_t
startTime_nEqBarConstr
=
time
(
NULL
);
localSum
=
0
.
0
;
long
off3
;
long
off2
;
for
(
int
ix
=
0
;
ix
<
nEqBarConstr
;
ix
++
){
//stars
for
(
int
ix
=
0
;
ix
<
nEqBarConstr
;
ix
++
)
{
//stars
yi
=
knownTerms
[
mapNoss
[
myid
]
+
nEqExtConstr
+
ix
];
long
offBarStarConstrEq
=
mapNcoeff
[
myid
]
+
nEqExtConstr
*
nOfElextObs
+
ix
*
nOfElBarObs
;
//Star offset on the row i2 (all the row)
#pragma omp for
for
(
long
yx
=
0
;
yx
<
numOfBarStar
;
yx
++
){
//FV_ EDIT ompSs
//#pragma omp for
for
(
long
yx
=
0
;
yx
<
numOfBarStar
;
yx
++
)
{
off3
=
yx
*
nAstroPSolved
;
off2
=
offBarStarConstrEq
+
off3
;
for
(
int
j2
=
0
;
j2
<
nAstroPSolved
;
j2
++
){
for
(
int
j2
=
0
;
j2
<
nAstroPSolved
;
j2
++
)
{
localSum
=
systemMatrix
[
off2
+
j2
]
*
yi
;
vVect
[
j2
+
off3
]
+=
localSum
;
}
}
#pragma omp barrier
//FV_ EDIT ompSs
//#pragma omp barrier
}
//for i2
if
(
comlsqr
.
itn
<=
2
&&
(
myid
==
0
||
myid
==
nproc
-
1
||
debugMode
==
1
))
printf
(
"
\t
PE=%d AprodTiming: nEqBarConstr sec=%ld
\n
"
,
myid
,
time
(
NULL
)
-
startTime_nEqBarConstr
);
}
//////////////////////////////////////////////////////
/// Mode 2 InstrConstr
if
(
nOfInstrConstr
){
if
(
nOfInstrConstr
)
{
time_t
startTime_nOfInstrConstr
=
time
(
NULL
);
localSum
=
0
.
0
;
int
offSetInstr
;
long
off1
;
...
...
@@ -433,7 +698,8 @@ if(nOfInstrConstr){
long
off4
=
mapNoss
[
myid
]
*
nInstrPSolved
;
long
off5
;
long
off6
;
for
(
int
k1
=
myid
;
k1
<
nOfInstrConstr
;
k1
=
k1
+
nproc
){
for
(
int
k1
=
myid
;
k1
<
nOfInstrConstr
;
k1
=
k1
+
nproc
)
{
yi
=
knownTerms
[
off2
+
k1
];
offSetInstr
=
0
;
for
(
int
m
=
0
;
m
<
k1
;
m
++
)
...
...
@@ -441,19 +707,24 @@ if(nOfInstrConstr){
off1
=
off3
+
offSetInstr
;
off5
=
off4
+
offSetInstr
;
#pragma omp for
for
(
int
j
=
0
;
j
<
instrConstrIlung
[
k1
];
j
++
){
//FV_ EDIT ompSs
//#pragma omp for
for
(
int
j
=
0
;
j
<
instrConstrIlung
[
k1
];
j
++
)
{
localSum
=
systemMatrix
[
off1
+
j
]
*
yi
;
off6
=
offInstrUnk
+
instrCol
[
off5
+
j
];
vVect
[
off6
]
+=
localSum
;
}
}
if
(
comlsqr
.
itn
<=
2
&&
(
myid
==
0
||
myid
==
nproc
-
1
||
debugMode
==
1
))
printf
(
"
\t
PE=%d AprodTiming: nOfInstrConstr sec=%ld
\n
"
,
myid
,
time
(
NULL
)
-
startTime_nOfInstrConstr
);
}
////////////
}
//pragma
*
ompSec
+=
time
(
NULL
)
-
startTime
;
if
(
comlsqr
.
itn
<=
2
&&
(
myid
==
0
||
myid
==
nproc
-
1
||
debugMode
==
1
))
printf
(
"PE=%d AprodTiming: mode=2.2 OmpSec=%ld
\n
"
,
myid
,
time
(
NULL
)
-
startTime
);
if
(
comlsqr
.
itn
<=
2
&&
(
myid
==
0
||
myid
==
nproc
-
1
||
debugMode
==
1
))
printf
(
"PE=%d AprodTiming: mode=2.2 OmpSec=%ld
\n
"
,
myid
,
time
(
NULL
)
-
startTime
);
}
// else if(mode==2)
}
This diff is collapsed.
Click to expand it.
lsqr.c
+
1314
−
1197
View file @
378fbe40
This diff is collapsed.
Click to expand it.
solvergaiaSim.c
+
2575
−
2194
View file @
378fbe40
This diff is collapsed.
Click to expand it.
Preview
0%
Loading
Try again
or
attach a new file
.
Cancel
You are about to add
0
people
to the discussion. Proceed with caution.
Finish editing this message first!
Save comment
Cancel
Please
register
or
sign in
to comment