feat(sparse): faster AddKel for sparse matrix + indent

labmec · May 27, 2024 · 051d75f · 051d75f
1 parent 1fa2090
commit 051d75f
Showing 1 changed file with 71 additions and 68 deletions.
diff --git a/Matrix/TPZYSMPMatrix.cpp b/Matrix/TPZYSMPMatrix.cpp
@@ -189,88 +189,91 @@ void TPZFYsmpMatrix<TVar>::AddKel(TPZFMatrix<TVar> & elmat, TPZVec<int64_t> & de
 }
 
 template<class TVar>
-void TPZFYsmpMatrix<TVar>::AddKel(TPZFMatrix<TVar> & elmat, TPZVec<int64_t> & sourceindex, TPZVec<int64_t> & destinationindex){
-	int64_t i,j,k = 0;
-	TVar value=0.;
-	int64_t ipos,jpos;
-	for(i=0;i<sourceindex.NElements();i++){
-		for(j=0;j<sourceindex.NElements();j++){
-			ipos=destinationindex[i];
-			jpos=destinationindex[j];
-			value=elmat.GetVal(sourceindex[i],sourceindex[j]);
-            //cout << "j= " << j << endl;
+void TPZFYsmpMatrix<TVar>::AddKel(TPZFMatrix<TVar>&elmat, TPZVec<int64_t> &sourceindex,
+																	TPZVec<int64_t> &destinationindex){
 
-			//cout << "fIA[ipos] " << fIA[ipos] << "     fIA[ipos+1] " << fIA[ipos+1] << endl;
-			int flag = 0;
-			k++;
-			if(k >= fIA[ipos] && k < fIA[ipos+1] && fJA[k]==jpos)
-			{ // OK -> elements in sequence
-			  fA[k]+=value;
-			  flag = 1;
-			}else
-			{
-			  for(k=fIA[ipos];k<fIA[ipos+1];k++){
-				if(fJA[k]==jpos || fJA[k]==-1){
-				  //cout << "fJA[k] " << fJA[k] << " jpos "<< jpos << "   " << value << endl;
-				  //cout << "k " << k << "   "<< jpos << "   " << value << endl;
-				  flag=1;
-				  if(fJA[k]==-1){
-					fJA[k]=jpos;
-					fA[k]=value;
-					// cout << jpos << "   " << value << endl;
-					break;
-				  }else{
-					fA[k]+=value;
-					break;
-				  }
+  // initialize original index locations
+  const int64_t neq = destinationindex.size();
+  TPZManVector<int64_t,800> idx(neq);
+  std::iota(idx.begin(), idx.end(), 0);
+
+  // sort indexes based on comparing values in v
+  // using std::stable_sort instead of std::sort
+  // to avoid unnecessary index re-orderings
+  // when v contains elements of equal values 
+  std::sort(idx.begin(), idx.end(),
+						[&destinationindex](const int64_t i1, const int64_t i2)
+						{return destinationindex[i1] < destinationindex[i2];});
+
+
+  for(auto dummy_i=0;dummy_i<neq;dummy_i++){
+    const auto i = idx[dummy_i];
+    const auto ipos=destinationindex[i];
+    //first col of the line
+    auto k = fIA[ipos];
+    const auto maxj = fIA[ipos+1];
+    const auto source_i = sourceindex[i];
+    for(auto dummy_j=0;dummy_j<neq;dummy_j++){
+      const auto j = idx[dummy_j];
+      const auto jpos=destinationindex[j];
+      const auto source_j = sourceindex[j];
+      const auto &value=elmat.GetVal(source_i,source_j);
+      while(fJA[k]!=jpos){
+				k++;
+				if(k==maxj){
+					std::cout << "TPZFYsmpMatrix::AddKelAtomic: "
+										<<" Non existing position on sparse matrix: "
+										<<" line =" << ipos << " column =" << jpos << endl;        
+					DebugStop();
 				}
-			  }
 			}
-			if(!flag) cout << "TPZFYsmpMatrix::AddKel: Non existing position on sparse matrix: line =" << ipos << " column =" << jpos << endl;         
-		}
-	}
+			fA[k]+=value;
+    }
+  }
 }
 
 template<class TVar>
 void TPZFYsmpMatrix<TVar>::AddKelAtomic(TPZFMatrix<TVar>&elmat, TPZVec<int64_t> &sourceindex,
-																				TPZVec<int64_t> &destinationindex){
+                                        TPZVec<int64_t> &destinationindex){
 
-	// initialize original index locations
-	const int64_t neq = destinationindex.size();
-	TPZManVector<int64_t,800> idx(neq);
-	std::iota(idx.begin(), idx.end(), 0);
+  // initialize original index locations
+  const int64_t neq = destinationindex.size();
+  TPZManVector<int64_t,800> idx(neq);
+  std::iota(idx.begin(), idx.end(), 0);
 
   // sort indexes based on comparing values in v
   // using std::stable_sort instead of std::sort
   // to avoid unnecessary index re-orderings
   // when v contains elements of equal values 
-	std::sort(idx.begin(), idx.end(),
-       [&destinationindex](size_t i1, size_t i2)
-			 {return destinationindex[i1] < destinationindex[i2];});
-
-
-	for(auto dummy_i=0;dummy_i<neq;dummy_i++){
-		const int64_t i = idx[dummy_i];
-		const auto ipos=destinationindex[i];
-		//first col of the line
-		int64_t k = fIA[ipos];
-		const auto maxj = fIA[ipos+1];
-		const auto source_i = sourceindex[i];
-		for(auto dummy_j=0;dummy_j<neq;dummy_j++){
-			const int j = idx[dummy_j];
-			const auto jpos=destinationindex[j];
-			const auto source_j = sourceindex[j];
-			const auto value=elmat.GetVal(source_i,source_j);
-			while(k<maxj && fJA[k]!=jpos){k++;}
-			if(k==maxj){
-				std::cout << "TPZFYsmpMatrix::AddKelAtomic: "
-									<<" Non existing position on sparse matrix: "
-									<<" line =" << ipos << " column =" << jpos << endl;        
-				DebugStop();
+  std::sort(idx.begin(), idx.end(),
+						[&destinationindex](const int64_t i1, const int64_t i2)
+						{return destinationindex[i1] < destinationindex[i2];});
+
+
+  for(auto dummy_i=0;dummy_i<neq;dummy_i++){
+    const auto i = idx[dummy_i];
+    const auto ipos=destinationindex[i];
+    //first col of the line
+    auto k = fIA[ipos];
+    const auto maxj = fIA[ipos+1];
+    const auto source_i = sourceindex[i];
+    for(auto dummy_j=0;dummy_j<neq;dummy_j++){
+      const auto j = idx[dummy_j];
+      const auto jpos=destinationindex[j];
+      const auto source_j = sourceindex[j];
+      const auto &value=elmat.GetVal(source_i,source_j);
+      while(fJA[k]!=jpos){
+				k++;
+				if(k==maxj){
+					std::cout << "TPZFYsmpMatrix::AddKelAtomic: "
+										<<" Non existing position on sparse matrix: "
+										<<" line =" << ipos << " column =" << jpos << endl;        
+					DebugStop();
+				}
 			}
-			pzutils::AtomicAdd(fA[k],value);
-		}
-	}
+      pzutils::AtomicAdd(fA[k],value);
+    }
+  }
 }
 
 template<class TVar>