# Code optimization

**URL:** <https://chapel.discourse.group/t/code-optimization/17990>\
**Category:** Users\
**Created:** [December 8, 2022, 10:12pm UTC](https://chapel.discourse.group/t/code-optimization/17990 "2022-12-08T22:12:24Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![annacfs](https://yyz2.discourse-cdn.com/free1/user_avatar/chapel.discourse.group/annacfs/32/361_2.png) [@annacfs](https://chapel.discourse.group/u/annacfs)\
**Post date:** [December 8, 2022, 10:12pm UTC](https://chapel.discourse.group/t/code-optimization/17990/1 "2022-12-08T22:12:24Z")

</div>

Hello guys,

I'm writing code to solve the navier stokes equations but the poisson problem isn't scaled as I'd like. Can anyone help me optimize the code?

Attached is a brief description and the code: [PDF description](https://drive.google.com/file/d/1aR0D4u72dYoPspcvSJNwLfX3jPe0Ta_n/view?usp=sharing)

> **[Google Drive: Sign-in](https://accounts.google.com/v3/signin/identifier?dsh=S-1609646196%3A1670537484315919&continue=https%3A%2F%2Fdrive.google.com%2Fdrive%2Ffolders%2F1FaUCXfWGBPrl-lgnlyGVGyjhHWy0YSuD%3Fusp%3Dsharing&followup=https%3A%2F%2Fdrive.google.com%2Fdrive%2Ffolders%2F1FaUCXfWGBPrl-lgnlyGVGyjhHWy0YSuD%3Fusp%3Dsharing&osid=1&passive=1209600&service=wise&flowName=GlifWebSignIn&flowEntry=ServiceLogin&ifkv=AeAAQh7eskblIW23LRvIACA3zLSXp61TIBccSG5eiREIEDUKsrBIuxatHQ5G4jT93df_QQ77BCjA5g)**
>
> Access Google Drive with a Google account (for personal use) or Google Workspace account (for business use).

I'm using single locale;

CHPL\_HOST\_PLATFORM: linux64  
CHPL\_HOST\_COMPILER: gnu  
CHPL\_HOST\_CC: gcc  
CHPL\_HOST\_CXX: g++  
CHPL\_HOST\_ARCH: x86\_64  
CHPL\_TARGET\_PLATFORM: linux64  
CHPL\_TARGET\_COMPILER: gnu  
CHPL\_TARGET\_CC: gcc  
CHPL\_TARGET\_CXX: g++  
CHPL\_TARGET\_ARCH: x86\_64  
CHPL\_TARGET\_CPU: native  
CHPL\_LOCALE\_MODEL: flat  
CHPL\_GPU\_CODEGEN: none  
CHPL\_COMM: none  
CHPL\_TASKS: fifo \*  
CHPL\_LAUNCHER: none  
CHPL\_TIMERS: generic  
CHPL\_UNWIND: none  
CHPL\_HOST\_MEM: cstdlib  
CHPL\_MEM: cstdlib \*  
CHPL\_ATOMICS: cstdlib  
CHPL\_GMP: none  
CHPL\_HWLOC: none \*  
CHPL\_RE2: none  
CHPL\_LLVM: unset  
CHPL\_LLVM\_CONFIG: none  
CHPL\_LLVM\_VERSION: 0  
CHPL\_AUX\_FILESYS: none  
CHPL\_LIB\_PIC: none  
CHPL\_SANITIZE: none  
CHPL\_SANITIZE\_EXE: none

module load chapel/1.27.0

source /modules/apps/chapel/chapel-1.27.0/util/setchplenv.bash && export CHPL\_MEM=cstdlib CHPL\_HWLOC=none CHPL\_TASKS=fifo CHPL\_RT\_NUM\_THREADS\_PER\_LOCALE=2

cd /lustre/anna2/poisson/chapel/grid130/

chpl --fast eliptico3D\_coforall.chpl

time ./eliptico3D\_coforall --imax=130 --jmax=130 --kmax=130 --npx=1 --npy=1 --npz=2 \> coforall\_130\_112

chpl --fast eliptico3D\_forall.chpl

time ./eliptico3D\_forall --imax=130 --jmax=130 --kmax=130 \> forall\_130\_2

Serial

```Chapel
// Finite Difference Solver for the Poisson equation (3D)  
// Author: JESUS, Anna && GRION, Livia   
// Last Modification: 20/11/2022
// Serial version
 
	use Time;
	var runtime: Timer;
	
	config const imax, jmax, kmax : int;	 
	const Lx = 1.0, Ly = 1.0, Lz = 1.0 : real;
	const dx = Lx/imax, dy = Ly/jmax, dz = Lz/kmax : real;      
	const w = 1.3 : real;										
	const dx2 = dx*dx, dy2 = dy*dy, dz2 = dz*dz: real;
	const beta1 = dx2/dy2 : real;	
	const beta2 = dx2/dz2 : real;	 
	const tol = 1E-6: real;
    const NPRIN = 500: int;
    
	var f, tf, residuo : [1..imax, 1..jmax, 1..kmax] real; 	  
	var i, j, k: int;
	var itera: int;
    var x, y, z, erro : real;
    var max_f, min_f : real;
    
 	runtime.start(); 
 	
 	serial{
    itera = 0;
    residuo = 0.0; 
    

	  for k in 1..kmax do{    
	   z = (k-0.5)*dz ;          
		
		for j in 1..jmax do{
		  y = (j-0.5)*dy;
		  
			for i in 1..imax do{ 
			    x = (i-0.5)*dx;	       		           		    
     		    
                tf(i,j,k) = -3.0*sin(x)*cos(y)*sin(z);
     		    
     		    }
     		}
        } 
        
            // Direcao - x
     		for j in 1..jmax do{     		   
     		    y = (j-0.5)*dy;       		    
     		        for k in 1..kmax do{
     		        z = (k-0.5)*dz;     		    
     		       
     		        x = 0.5;
                    f(1,j,k) = sin(x)*cos(y)*sin(z);
                
                    x = (imax-0.5)*dx;
                    f(imax,j,k) = sin(x)*cos(y)*sin(z);
     		         
     		    }
     	    }
     		    
     		     
            // Direcao - y 
	        for i in 1..imax do{ 
	            x = (i-0.5)*dx;
     		    
     		    for k in 1..kmax do{
     		        z = (k-0.5)*dy;
     		    
     		        y = 0.5;
                    f(i,1,k) = sin(x)*cos(y)*sin(z); 
                    
                    y = (jmax-0.5)*dy;
                    f(i,jmax,k) = sin(x)*cos(y)*sin(z);  
     		     
     		    }
     		 }
     		     
         // Direcao - z 	 	   
	     for i in 1..imax do{
	        x = (i-0.5)*dx; 	       
     		   
     		   for j in 1..jmax do{ 
     		   y = (j-0.5)*dy; 
     		    
     		   z = 0.5;          
               f(i,j,1) = sin(x)*cos(y)*sin(z);
                
               z = (kmax-0.5)*dz ;        
               f(i,j,kmax) = sin(x)*cos(y)*sin(z);  
     		    
     		     
     		 }
        }
     		     
		erro = 1E0;
		itera = 0;
  
		    while (erro >= tol) do{    
            
		    f = sor_method(f, tf, imax, jmax, kmax, dx2, beta1, beta2, w);		 			 		 
			residuo = calc_erro(f, tf, imax, jmax, kmax, dx2, dy2, dz2);   

            erro = max reduce(residuo); 

            itera += 1;	  
					
		    
		        if (mod(itera,NPRIN) == 0) then{		  		
		            writef('%i %6.12dr\n', itera, erro); 
		        } 
  
            }
        } // fim serial
 		 
 		runtime.stop();
 		 	
   		escreve(f, imax, jmax, kmax, dx, dy, dz);
   		
 
        max_f = max reduce(f); 
        min_f = min reduce(f);     

	writef("%12.8er %12.8er \n", max_f, min_f);				                 				     
    writef(" CPU TIME CHAPEL = %10.16dr \n",runtime.elapsed());
	writef("%i %12.8er \n",itera, erro);						                 
	writef("------------------------------------------------------------------------------ \n");  
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////
  
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////
	proc sor_method(f, tf, imax: int, jmax: int, kmax: int, dx2: real, beta1: real, beta2: real, w: real){  
	
     for i in 2..imax-1 do{
     		for j in 2..jmax-1 do{
     		    for k in 2..kmax-1 do{
     		  
                    f(i,j,k) = (1.0 - w)*f(i,j,k) +  
                                w * (-dx2*tf(i,j,k) + f(i-1,j,k) + f(i+1,j,k) +
                                beta1*( f(i,j-1,k) + f(i,j+1,k) ) + 
                                beta2*( f(i,j,k-1) + f(i,j,k+1) ))/ ( 2.0*(beta1 + beta2 + 1.0) ) ;
            }   
        }   
    }   
 
		return f;     
       
}
////////////////////////////////////////////////////////////////////////////////////////////////////////////////// 

	proc calc_erro(f, tf, imax: int, jmax: int, kmax: int, dx2: real, dy2: real, dz2: real){

	var errol : [1..imax,1..jmax, 1..kmax] real; 
	errol = 0.0;

	 for i in 2..imax-1 do{
     		for j in 2..jmax-1 do{
     		    for k in 2..kmax-1 do{
		 
      		        errol(i,j,k) = abs(tf(i,j,k) -  
			                               (f(i-1,j,k) - 2.0*f(i,j,k) + f(i+1,j,k))/dx2 -
			                               (f(i,j-1,k) - 2.0*f(i,j,k) + f(i,j+1,k))/dy2 -
			                               (f(i,j,k-1) - 2.0*f(i,j,k) + f(i,j,k+1))/dz2);   
 		    }
  	    }
  	}
     		return errol;

}

 ////////////////////////////////////////////////////////////////////////////////////////////////////////////////// 
    proc escreve(f, imax: int, jmax: int, kmax: int, dx: real, dy: real, dz: real){  
 
    use IO;                  
	var erro: [1..imax,1..jmax, 1..kmax] real; 
	erro = 0.0;
    const saida = openwriter("saida_serial.csv");  
   
    saida.writef(" %s %s %s %s %s %s %s %s %s \n", "x", ",", "y", ",", "z", ",", "f", ",", "e");
	for k in 1..kmax do{    
	   z = (k-0.5)*dz ;          
		
		for j in 1..jmax do{
		  y = (j-0.5)*dy;
		  
			for i in 1..imax do{ 
			    x = (i-0.5)*dx;
			            
			    erro(i,j,k) = abs(f(i,j,k) - sin(x)*cos(y)*sin(z));
			            
        		saida.writef("%10.12dr %s %10.12dr %s %10.12dr %s %10.12dr %s %10.12dr \n", 
        		            x, ",", y, ",", z, ",", f(i,j,k), ",", erro(i,j,k));
        	 }
 
		  }
    }
saida.close();
 
 }
 

```

Coforall

```Chapel
// Finite Difference Solver for the Poisson equation (3D)  
// Author: JESUS, Anna && GRION, Livia   
// Last Modification: 20/11/2022
// Serial version
 
	use Time;
	var runtime: Timer;
	
	config const imax, jmax, kmax : int;	 
    config const npx, npy, npz: int;	  
    const ptsx = (imax-2)/npx, ptsy = (jmax-2)/npy, ptsz = (kmax-2)/npz : int;    
	const Lx = 1.0, Ly = 1.0, Lz = 1.0 : real;
	const dx = Lx/imax, dy = Ly/jmax, dz = Lz/kmax : real;      
	const w = 1.3 : real;										
	const dx2 = dx*dx, dy2 = dy*dy, dz2 = dz*dz: real;
	const beta1 = dx2/dy2 : real;	
	const beta2 = dx2/dz2 : real;	 
	const tol = 1E-6: real;
    const NPRIN = 500;
    const Inner = {2..imax-1, 2..jmax-1, 2..kmax-1};
	
	var f, tf, residuo : [1..imax, 1..jmax, 1..kmax] real; 	  
 
	var i, j, k: int;
	var itera: int;
    var x, y, z, erro : real;
    var max_f, min_f : real;
    
 	runtime.start(); 
 	
    itera = 0;
    residuo = 0.0;
 
	  for k in 1..kmax do{    
	   z = (k-0.5)*dz ;          
		
		for j in 1..jmax do{
		  y = (j-0.5)*dy;
		  
			for i in 1..imax do{ 
			    x = (i-0.5)*dx;
     		    
                tf(i,j,k) = -3.0*sin(x)*cos(y)*sin(z);
     		    
     		    }
     		}
       }
        
            // Direcao - x
     		for j in 1..jmax do{     		   
     		    y = (j-0.5)*dy;       		    
     		        for k in 1..kmax do{
     		        z = (k-0.5)*dz;     		    
     		       
     		        x = 0.5;
                    f(1,j,k) = sin(x)*cos(y)*sin(z);
                
                    x = (imax-0.5)*dx;
                    f(imax,j,k) = sin(x)*cos(y)*sin(z);
     		         
     		    }
     	    }
     		     
            // Direcao - y 
	        for i in 1..imax do{ 
	            x = (i-0.5)*dx;
     		    
     		    for k in 1..kmax do{
     		        z = (k-0.5)*dy;
     		    
     		        y = 0.5;
                    f(i,1,k) = sin(x)*cos(y)*sin(z); 
                    
                    y = (jmax-0.5)*dy;
                    f(i,jmax,k) = sin(x)*cos(y)*sin(z);  
     		     
     		    }
     		 }
     		     
         // Direcao - z  	   
	     for i in 1..imax do{
	        x = (i-0.5)*dx; 	       
     		   
     		   for j in 1..jmax do{ 
     		   y = (j-0.5)*dy; 
     		    
     		   z = 0.5;          
               f(i,j,1) = sin(x)*cos(y)*sin(z);
                
               z = (kmax-0.5)*dz ;        
               f(i,j,kmax) = sin(x)*cos(y)*sin(z);  
     		    
     		     
     		 }
        }
	 
		erro = 1E0;
		itera = 0;
 
		    while (erro >= tol) do{    
            
		    f = sor_method(f, tf, imax, jmax, kmax, dx2, beta1, beta2, w, ptsx, ptsy, ptsz);		 			 		 
			residuo = calc_erro(f, tf, imax, jmax, kmax, dx2, dy2, dz2, ptsx, ptsy, ptsz);   

            erro = max reduce(residuo); 

            itera += 1;	  
							    
		   if (mod(itera,NPRIN) == 0) then{		  		
		    writef('%i %6.12dr\n', itera, erro); 
		   } 
  
        }
 		 
 		runtime.stop();
 		 	
   	    escreve(f, imax, jmax, kmax, dx, dy, dz); 
    
    max_f = max reduce(f); 
    min_f = min reduce(f); 
    
    writef("%12.8er %12.8er \n", max_f, min_f);				                 				     
    writef(" CPU TIME CHAPEL = %10.16dr \n",runtime.elapsed());
	writef("%i %12.8er \n",itera, erro);						                 
	writef("------------------------------------------------------------------------------ \n");  
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////
  
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////
	proc sor_method(f, tf, imax: int, jmax: int, kmax: int, dx2: real, beta1: real, beta2: real, w: real, ptsx: int, ptsy: int, ptsz: int){ 
  
    
 	coforall (ii,jj,kk) in {2..imax-1, 2..jmax-1, 2..kmax-1} by (ptsx, ptsy, ptsz) {
 
     for i in ii..ii + ptsx-1 do{
     		for j in jj..jj + ptsy-1 do{
     		    for k in kk..kk + ptsz-1 do{
     		  
                    f(i,j,k) = (1.0 - w)*f(i,j,k) +  
                                w * (-dx2*tf(i,j,k) + f(i-1,j,k) + f(i+1,j,k) +
                                beta1*( f(i,j-1,k) + f(i,j+1,k) ) + 
                                beta2*( f(i,j,k-1) + f(i,j,k+1) ))/ ( 2.0*(beta1 + beta2 + 1.0) ) ;
            }   
        }   
    } 
    }  
 
		return f; 
       
}
////////////////////////////////////////////////////////////////////////////////////////////////////////////////// 

	proc calc_erro(f, tf, imax: int, jmax: int, kmax: int, dx2: real, dy2: real, dz2: real, ptsx: int, ptsy: int, ptsz: int){

	var errol : [1..imax,1..jmax, 1..kmax] real; 
	errol = 0.0;

	coforall (ii,jj,kk) in {2..imax-1, 2..jmax-1, 2..kmax-1} by (ptsx, ptsy, ptsz) {
 
     for i in ii..ii + ptsx-1 do{
     		for j in jj..jj + ptsy-1 do{
     		    for k in kk..kk + ptsz-1 do{
     		  
		 
      		        errol(i,j,k) = abs(tf(i,j,k) -  
			                               (f(i-1,j,k) - 2.0*f(i,j,k) + f(i+1,j,k))/dx2 -
			                               (f(i,j-1,k) - 2.0*f(i,j,k) + f(i,j+1,k))/dy2 -
			                               (f(i,j,k-1) - 2.0*f(i,j,k) + f(i,j,k+1))/dz2);   
 		    }
  	    }
  	}
     
     }
		return errol;

}
 

 ////////////////////////////////////////////////////////////////////////////////////////////////////////////////// 
    proc escreve(f, imax: int, jmax: int, kmax: int, dx: real, dy: real, dz: real){  
 
    use IO;                  
	var erro: [1..imax,1..jmax, 1..kmax] real; 
	erro = 0.0; 
	
    const saida = openwriter("saida_coforall.csv");  
   
    saida.writef(" %s %s %s %s %s %s %s %s %s \n", "x", ",", "y", ",", "z", ",", "f", ",", "e");
	for k in 1..kmax do{    
	   z = (k-0.5)*dz ;          
		
		for j in 1..jmax do{
		  y = (j-0.5)*dy;
		  
			for i in 1..imax do{ 
			    x = (i-0.5)*dx;
			            
			    erro(i,j,k) = abs(f(i,j,k) - sin(x)*cos(y)*sin(z));
			            
        		saida.writef("%10.12dr %s %10.12dr %s %10.12dr %s %10.12dr %s %10.12dr \n", 
        		            x, ",", y, ",", z, ",", f(i,j,k), ",", erro(i,j,k));
        	 }
 
		  }
    }
saida.close();
 
 }

```

Forall

```Chapel
// Finite Difference Solver for the Poisson equation (3D)  
// Author: JESUS, Anna && GRION, Livia   
// Last Modification: 20/10/2022
// Parallel version
 
	use Time;
	var runtime: Timer;
	
	config const imax, jmax, kmax : int;	 
	const Lx = 1.0, Ly = 1.0, Lz = 1.0 : real;
	const dx = Lx/imax, dy = Ly/jmax, dz = Lz/kmax : real;       
	const w = 1.3 : real;										
	const dx2 = dx*dx, dy2 = dy*dy, dz2 = dz*dz: real;
	const beta1 = dx2/dy2 : real;	
	const beta2 = dx2/dz2 : real;	 
	const tol = 1E-6: real;
 	const Inner = {2..imax-1, 2..jmax-1, 2..kmax-1};
    const domPoisson = {1..imax, 1..jmax, 1..kmax};
    const NPRIN = 500;
     		
	var f, tf, residuo : [1..imax, 1..jmax, 1..kmax] real; 	  
	var i, j, k: int;
	var itera: int;
    var x, y, z, erro : real; 
    var max_f, min_f : real;
        
        
    itera = 0;
    residuo = 0.0;
	
	runtime.start();
	 
	   for k in 1..kmax do{    
	   z = (k-0.5)*dz ;          
		
		for j in 1..jmax do{
		  y = (j-0.5)*dy;
		  
			for i in 1..imax do{ 
			    x = (i-0.5)*dx;  		     
	       		    
                tf(i,j,k) = -3.0*sin(x)*cos(y)*sin(z);      		    
            
            }
          }
        
        }
        
            // Direcao - x
     	    for j in 1..jmax do{     		   
     		    y = (j-0.5)*dy;       		    
     		        for k in 1..kmax do{
     		        z = (k-0.5)*dz;     		    
     		       
     		        x = 0.5;
                    f(1,j,k) = sin(x)*cos(y)*sin(z);
                
                    x = (imax-0.5)*dx;
                    f(imax,j,k) = sin(x)*cos(y)*sin(z);
     		         
     		    }
     	    } 
     		     
            // Direcao - y 
	        for i in 1..imax do{ 
	            x = (i-0.5)*dx;
     		    
     		    for k in 1..kmax do{
     		        z = (k-0.5)*dy;
     		    
     		        y = 0.5;
                    f(i,1,k) = sin(x)*cos(y)*sin(z); 
                    
                    y = (jmax-0.5)*dy;
                    f(i,jmax,k) = sin(x)*cos(y)*sin(z);  
     		     
     		    }
     		 }
     		      
     		   
	     for i in 1..imax do{
	        x = (i-0.5)*dx; 	       
     		   
     		   for j in 1..jmax do{ 
     		   y = (j-0.5)*dy; 
     		    
     		   z = 0.5;          
               f(i,j,1) = sin(x)*cos(y)*sin(z);
                
               z = (kmax-0.5)*dz ;        
               f(i,j,kmax) = sin(x)*cos(y)*sin(z);  
     		    
     		     
     		 }
        } 
        
		erro = 1E0;
		itera = 0;
  
		    while (erro >= tol) do{    
            
		    f = sor_method(f, tf, dx2, beta1, beta2, w);		 		 
			residuo = calc_erro(f, tf, dx2, dy2, dz2);   
			erro = max reduce(residuo); 
            itera += 1;	  			
		  
		        if (mod(itera,NPRIN) == 0) then{		  		
		            writef('%i %6.12dr\n', itera, erro); 
		        } 
        
           } 
		  
   		 escreve(f, imax, jmax, kmax, dx, dy, dz);
  
	runtime.stop(); 
	
    max_f = max reduce(f); 
    min_f = min reduce(f);  
	 	
	writef("%12.8er %12.8er \n", max_f, min_f);				                 				     
    writef(" CPU TIME CHAPEL = %10.16dr \n",runtime.elapsed());
	writef("%i %12.8er \n",itera, erro);						                 
	writef("------------------------------------------------------------------------------ \n"); 
////////////////////////////////////////////////////////////////////////////////////////////////////////////////// 

//////////////////////////////////////////////////////////////////////////////////////////////////////////////////
	proc sor_method(f, tf, dx2: real, beta1: real, beta2: real, w: real){ 
 
 
     forall (i,j,k) in Inner do{
     		  
                    f(i,j,k) = (1.0 - w)*f(i,j,k) +  
                                w * (-dx2*tf(i,j,k) + f(i-1,j,k) + f(i+1,j,k) +
                                beta1*( f(i,j-1,k) + f(i,j+1,k) ) + 
                                beta2*( f(i,j,k-1) + f(i,j,k+1) ))/ ( 2.0*(beta1 + beta2 + 1.0) ) ;
            
    }    
		return f;     
       
} 
////////////////////////////////////////////////////////////////////////////////////////////////////////////////// 

	proc calc_erro(f, tf, dx2: real, dy2: real, dz2: real){

	var errol : [1..imax,1..jmax, 1..kmax] real; 
	errol = 0.0;

	forall (i,j,k) in Inner do{ 
	
      		        errol(i,j,k) = abs(tf(i,j,k) -  
			                               (f(i-1,j,k) - 2.0*f(i,j,k) + f(i+1,j,k))/dx2 -
			                               (f(i,j-1,k) - 2.0*f(i,j,k) + f(i,j+1,k))/dy2 -
			                               (f(i,j,k-1) - 2.0*f(i,j,k) + f(i,j,k+1))/dz2);    
  	    }
     
		return errol;

} 
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////

    proc escreve(f, imax: int, jmax: int, kmax: int, dx: real, dy: real, dz: real){  
 
    use IO;                  
	var erro: [1..imax,1..jmax, 1..kmax] real; 
    const saida = openwriter("saida_chapel_par.csv");  
   
   
    saida.writef(" %s %s %s %s %s %s %s %s %s \n", "x", ",", "y", ",", "z", ",", "f", ",", "e");
	for k in 1..kmax do{    
	   z = (k-0.5)*dz ;          
		
		for j in 1..jmax do{
		  y = (j-0.5)*dy;
		  
			for i in 1..imax do{ 
			    x = (i-0.5)*dx;
			            
			    erro(i,j,k) = abs(f(i,j,k) - sin(x)*cos(y)*sin(z));
			            
        		saida.writef("%10.12dr %s %10.12dr %s %10.12dr %s %10.12dr %s %10.12dr \n", 
        		            x, ",", y, ",", z, ",", f(i,j,k), ",", erro(i,j,k));
        	 }
 
		  }
    }
saida.close();
 
 } 

```

---

<div class="post-metadata">

**Author:** ![damianmoz](https://yyz2.discourse-cdn.com/free1/user_avatar/chapel.discourse.group/damianmoz/32/319_2.png) [@damianmoz](https://chapel.discourse.group/u/damianmoz)\
**Post date:** [December 8, 2022, 10:23pm UTC](https://chapel.discourse.group/t/code-optimization/17990/2 "2022-12-08T22:23:44Z")

</div>

Anna, can you please attach the PDF and the raw chapel directly rather than require people to go to Google drive which requires separate access? Just a thought.

Did you implement all those earlier suggestions which were provided?

---

<div class="post-metadata">

**Author:** ![annacfs](https://yyz2.discourse-cdn.com/free1/user_avatar/chapel.discourse.group/annacfs/32/361_2.png) [@annacfs](https://chapel.discourse.group/u/annacfs)\
**Post date:** [December 8, 2022, 10:52pm UTC](https://chapel.discourse.group/t/code-optimization/17990/4 "2022-12-08T22:52:51Z")

</div>

what suggestions? I can not upload the files here

---

<div class="post-metadata">

**Author:** ![damianmoz](https://yyz2.discourse-cdn.com/free1/user_avatar/chapel.discourse.group/damianmoz/32/319_2.png) [@damianmoz](https://chapel.discourse.group/u/damianmoz)\
**Post date:** [December 8, 2022, 10:56pm UTC](https://chapel.discourse.group/t/code-optimization/17990/5 "2022-12-08T22:56:29Z")

</div>

OK. I did some cut and paste from the post and I now have the Chapel. I will see if I can compile it.

I will see what the Google drive request does..

---

<div class="post-metadata">

**Author:** ![damianmoz](https://yyz2.discourse-cdn.com/free1/user_avatar/chapel.discourse.group/damianmoz/32/319_2.png) [@damianmoz](https://chapel.discourse.group/u/damianmoz)\
**Post date:** [December 8, 2022, 11:00pm UTC](https://chapel.discourse.group/t/code-optimization/17990/6 "2022-12-08T23:00:27Z")

</div>

I grabbed the PDF also.

---

<div class="post-metadata">

**Author:** ![jcorrado](https://avatars.discourse-cdn.com/v4/letter/j/54ee81/32.png) [@jcorrado](https://chapel.discourse.group/u/jcorrado)\
**Post date:** [December 8, 2022, 11:28pm UTC](https://chapel.discourse.group/t/code-optimization/17990/7 "2022-12-08T23:28:32Z")

</div>

Hello @annacfs, thanks for the question!

One very quick change that should make a reasonable performance difference: remove the return statements from `sor_method` and `calc_erro`.

- In `sor_method`, `f` doesn't have to be returned because it's already an argument. (The array `f` will be modified by reference, so the `f = ` part of the call is also not needed. It can just be called as: `sor_method(f, tf, dx2, beta1, beta2, w);`)

- For `calc_error`, I'd delete the `errol` declaration and add an argument with the same name instead. So your proc would look like:

---

<div class="post-metadata">

**Author:** ![damianmoz](https://yyz2.discourse-cdn.com/free1/user_avatar/chapel.discourse.group/damianmoz/32/319_2.png) [@damianmoz](https://chapel.discourse.group/u/damianmoz)\
**Post date:** [December 9, 2022, 12:08am UTC](https://chapel.discourse.group/t/code-optimization/17990/8 "2022-12-09T00:08:24Z")

</div>

There were some suggestions made to the posts you made earlier in the year. I have not yet had time to look at the latest code in depth but it looks like you are still to implement all of those previous ideas.

Off the cuff,  
a) there are far too many divisions and reliance on the compiler to optimize expressions; and  
b) there are cases where a \*_forall_ should wrap the innermost loop which should use a **for**

I could be wrong but I think that you over parallelised the **sor\_method** which is an example of (b).

```chapel
forall (i, j) in do
{
    foreach k in 2..kmax-1 do
    {
       ...
    }
}

```

might be a better approach. I try to avoid having a **forall** as the innermost loop.

@jcorrado is on the right path - some of what he said was mentioned earlier in the year. The parameter t should/must have an **inout** intent. The overhead of making a copy would be horrendous.

Can you provide more detailed timings of individual segments of the program?

Thanks for the raw email.

---

<div class="post-metadata">

**Author:** ![damianmoz](https://yyz2.discourse-cdn.com/free1/user_avatar/chapel.discourse.group/damianmoz/32/319_2.png) [@damianmoz](https://chapel.discourse.group/u/damianmoz)\
**Post date:** [December 9, 2022, 12:23am UTC](https://chapel.discourse.group/t/code-optimization/17990/9 "2022-12-09T00:23:58Z")

</div>

Is Eq(2) recursive?

It seems that once you calculate f(i, j, 10) say, you need to use that value to then calculate f(i,j,11). Is that correct?

---

<div class="post-metadata">

**Author:** ![damianmoz](https://yyz2.discourse-cdn.com/free1/user_avatar/chapel.discourse.group/damianmoz/32/319_2.png) [@damianmoz](https://chapel.discourse.group/u/damianmoz)\
**Post date:** [December 9, 2022, 8:31am UTC](https://chapel.discourse.group/t/code-optimization/17990/10 "2022-12-09T08:31:29Z")

</div>

Consider the following. You might need to experiment. If `imax` is small, you might need to change this.

```Chapel
    proc sor_method(f, tf, dx2: real, beta1: real, beta2: real, w: real)
    {
        // recompute what is stored inside Inner to avoid
        // using data outside the scope of this routine
        const ijk = f.domain.dims();
        const _i = ijk[0].expand(-1); // 2..imax-1
        const _j = ijk[1].expand(-1); // 2..jmax-1
        const _k = ijk[2].expand(-1); // 2..kmax-1
        const w1 = (1.0 - w);
        const w2 = (0.5 * w) / (beta1 + beta2 + 1.0);

        forall i in _i do // if imax is small, this might be sub-optimal
        {
            // the entire loop contents will be run as a single task
            // so ensure that each such task has lots of work to do

            for j in _j do
            {
                const ref fijm1 = f[i, j-1, ..], fijp1 = f[i, j+1, ..];
                const ref fim1j = f[i-1, j, ..], fip1j = f[i+1, j, ..];
                const ref tfij = tf[i, j, ..];
                ref fij = f[i, j, ..];
                var fijk = f[i, j, 1];

                for k in _k do // recursive - cannot be vectorized - sigh!!
                {
                    // superscalar and FMA opportunities here
                    // this ordering might be suboptimal but

                    const b0 = fim1j[k] + fip1j[k];
                    const b1 = fijm1[k] + fijp1[k];
                    const b2 = fijk + fij[k+1];
                    const bb = b1 * beta1 + b2 * beta2 + b0 - dx2 * tfij[k];

                    fijk = w1 * fij[k] + w2 * bb;
                    fij[k] = fijk;
                }
            }
        }
    }

```

Use the experience from implementing this with @jcorrado's advice to create a better calc\_error

---

<div class="post-metadata">

**Author:** ![damianmoz](https://yyz2.discourse-cdn.com/free1/user_avatar/chapel.discourse.group/damianmoz/32/319_2.png) [@damianmoz](https://chapel.discourse.group/u/damianmoz)\
**Post date:** [December 9, 2022, 8:53am UTC](https://chapel.discourse.group/t/code-optimization/17990/11 "2022-12-09T08:53:30Z")

</div>

If vectorization is actually working well, the following might be better

```Chapel
    proc sor_kernel(_k, f, tf, dx2: real, beta1: real)
    {
        var b : [_k] real;
        const ref fijm1 = f[i, j-1, ..];
        const ref fijp1 = f[i, j+1, ..];
        const ref fim1j = f[i-1, j, ..];
        const ref fip1j = f[i+1, j, ..];
        const ref tfij = tf[i, j, ..];

        foreach k in _k do // this can be vectorized
        {
            // there are superscalar and FMA opportunities here
            // note this ordering might be suboptimal here goes

            const b1 = fijm1[k] + fijp1[k];
            const b0 = fim1j[k] + fip1j[k];

            b[k] = b1 * beta1 + b0 - dx2 * tfij[k];
        }
        return b;
    }
    proc sor_method(f, tf, dx2: real, beta1: real, beta2: real, w: real)
    {
        // recompute what is stored inside Inner to avoid
        // using data outside the scope of this routine
        const ijk = f.domain.dims();
        const _i = ijk[0].expand(-1); // 2..imax-1
        const _j = ijk[1].expand(-1); // 2..jmax-1
        const _k = ijk[2].expand(-1); // 2..kmax-1
        const w1 = (1.0 - w);
        const w2 = (0.5 * w) / (beta1 + beta2 + 1.0);

        forall i in _i do // if imax is small, this might be sub-optimal
        {
            // the entire loop contents will be run as a single task
            // so ensure that each such task has lots of work to do

            for j in _j do
            {
                ref fij = f[i, j, ..];
                var fijk = f[i, j, 1];
                const b = sor_kernel(_k, f, tf, dx2, beta1);

                for k in _k do // recursive - cannot be vectorized - sigh!!
                {
                    const b2 = fijk + fij[k+1];

                    fijk = w1 * fij[k] + w2 * (b2 * beta2 + b[k]);
                    fij[k] = fijk;
                }
            }
        }
    }

```

Try it. Experiment.

---

<div class="post-metadata">

**Author:** ![annacfs](https://yyz2.discourse-cdn.com/free1/user_avatar/chapel.discourse.group/annacfs/32/361_2.png) [@annacfs](https://chapel.discourse.group/u/annacfs)\
**Post date:** [December 9, 2022, 10:02am UTC](https://chapel.discourse.group/t/code-optimization/17990/12 "2022-12-09T10:02:56Z")

</div>

Nice!!!!! Let me think about it!  
Thank u

---

<div class="post-metadata">

**Author:** ![damianmoz](https://yyz2.discourse-cdn.com/free1/user_avatar/chapel.discourse.group/damianmoz/32/319_2.png) [@damianmoz](https://chapel.discourse.group/u/damianmoz)\
**Post date:** [December 9, 2022, 3:31pm UTC](https://chapel.discourse.group/t/code-optimization/17990/13 "2022-12-09T15:31:17Z")

</div>

I have not tested this. I hope I have not lost anything in the cutting and pasting I did. Others may have some better ideas. I am not sure how well the currrent version of the compiler is at vectorizing the code, something about which @mppf might have way better insight than I. Are you using the LLVM or GCC based variant of the compiler, and again, others can probably better advise you here but that will influence how well your code will vectorize. Maybe @jcorrado has some insightt here as well?

What are representative values of the problem size, i.e. what are the values of imax, jmax and kmax?

---

<div class="post-metadata">

**Author:** ![damianmoz](https://yyz2.discourse-cdn.com/free1/user_avatar/chapel.discourse.group/damianmoz/32/319_2.png) [@damianmoz](https://chapel.discourse.group/u/damianmoz)\
**Post date:** [December 9, 2022, 3:44pm UTC](https://chapel.discourse.group/t/code-optimization/17990/14 "2022-12-09T15:44:55Z")

</div>

See also the discussion on Github issue #19636 which is related.

---

<div class="post-metadata">

**Author:** ![annacfs](https://yyz2.discourse-cdn.com/free1/user_avatar/chapel.discourse.group/annacfs/32/361_2.png) [@annacfs](https://chapel.discourse.group/u/annacfs)\
**Post date:** [December 9, 2022, 6:03pm UTC](https://chapel.discourse.group/t/code-optimization/17990/15 "2022-12-09T18:03:39Z")

</div>

I tested imax=kmax=kmax = 130 in cluster. The size of the problem is on the order of 10^6; But my goal is to do 10^9

---

<div class="post-metadata">

**Author:** ![damianmoz](https://yyz2.discourse-cdn.com/free1/user_avatar/chapel.discourse.group/damianmoz/32/319_2.png) [@damianmoz](https://chapel.discourse.group/u/damianmoz)\
**Post date:** [December 9, 2022, 7:03pm UTC](https://chapel.discourse.group/t/code-optimization/17990/16 "2022-12-09T19:03:52Z")

</div>

Thanks. My idea of using the **forall** on the outside loop is OK then.

Note that if you ultimately want to run this on multiple nodes of a cluster, you might need to use **coforall** and tweak the design of **sor\_method**. But that can wait for now. Get things running nicely on your single node for now. I assume that the single node you are using in the cluster has two 28-core CPUs, meaning you have 56 cores available to you.

Also, in my opinion, your need to turn multi-threading off on your system. We always get better performance having a single thrread for each FPU. Others might have different experience using Chapel and I suggest you get their advice.

---

<div class="post-metadata">

**Author:** ![damianmoz](https://yyz2.discourse-cdn.com/free1/user_avatar/chapel.discourse.group/damianmoz/32/319_2.png) [@damianmoz](https://chapel.discourse.group/u/damianmoz)\
**Post date:** [December 9, 2022, 11:20pm UTC](https://chapel.discourse.group/t/code-optimization/17990/17 "2022-12-09T23:20:06Z")

</div>

You could/should also look at unrolling the innermost loops.

I would also be curious whether a **foreach** loop is faster than an **for** loop unrolled by 4. If the former is no faster, then the original suggestion, with its innermost loop unrolled by at least 2 (or maybe 4), would probably be optimal.

---

<div class="post-metadata">

**Author:** ![damianmoz](https://yyz2.discourse-cdn.com/free1/user_avatar/chapel.discourse.group/damianmoz/32/319_2.png) [@damianmoz](https://chapel.discourse.group/u/damianmoz)\
**Post date:** [December 10, 2022, 1:05am UTC](https://chapel.discourse.group/t/code-optimization/17990/18 "2022-12-10T01:05:25Z")

</div>

You should force imax, jmax,and kmax to be an odd multiple of 2, This means the number of elements in a range like

```Chapel
2..imax-1

```

is always a multiple of 4. This makes loop unrolling a lot easier to implement. My 2c.

---

<div class="post-metadata">

**Author:** ![damianmoz](https://yyz2.discourse-cdn.com/free1/user_avatar/chapel.discourse.group/damianmoz/32/319_2.png) [@damianmoz](https://chapel.discourse.group/u/damianmoz)\
**Post date:** [December 10, 2022, 1:57pm UTC](https://chapel.discourse.group/t/code-optimization/17990/19 "2022-12-10T13:57:18Z")

</div>

I ran your task on a 24-core system,. I expected a gain of about 11,

With the original code, I get  
a) a serial performance of 690 seconds in 15000+ iterations  
b) a parallel performance of 53 seconds in 22000+ iterations

The performance gain is about a factor of 13 which is slightly better that I thought. Do you numbers agree to within some tolerance? The parallelization of that recursive formula worries me.

When I run the modified version using the approach I suggested, I get  
b) a parallel performance of 73 seconds in 15000+ iterations

The performance gain is about 9.5 which is still pretty good and the iteration count matches.

EDITED **after the event**

I rewrote in a simpler style, still handling the recursive aspects of the 3rd dimension loop, and the performance gain is now nearly 22 on the 24 core machine (with a total of 8 memory channels) and nearly 28 for the 52-core machine (with a total of 12 memory channels). The figure of 28 is misleading because the serial version runs at 3.9Ghz if it is the only task on this CPU whereas the parallel version with 52 cores runs at 2.1Ghz. If you scale that up, you get a factor of about 50 on a 52 core system. I will post the code tomorrow.

That said, I think there is some issue the allocation of cores on your cluster.

---

<div class="post-metadata">

**Author:** ![damianmoz](https://yyz2.discourse-cdn.com/free1/user_avatar/chapel.discourse.group/damianmoz/32/319_2.png) [@damianmoz](https://chapel.discourse.group/u/damianmoz)\
**Post date:** [December 10, 2022, 3:04pm UTC](https://chapel.discourse.group/t/code-optimization/17990/20 "2022-12-10T15:04:01Z")

</div>

Can you please provide the Fortran+MPI reference code?

---

<div class="post-metadata">

**Author:** ![damianmoz](https://yyz2.discourse-cdn.com/free1/user_avatar/chapel.discourse.group/damianmoz/32/319_2.png) [@damianmoz](https://chapel.discourse.group/u/damianmoz)\
**Post date:** [December 10, 2022, 5:47pm UTC](https://chapel.discourse.group/t/code-optimization/17990/21 "2022-12-10T17:47:59Z")

</div>

What do your 2 graphs mean? What serial run time do you get?

When you say 28 cores, 56 threads, what are the actual CPUs, i.e.

```Chapel
cat /proc/cpuinfo

```

on the Linux machine. The machine I am using has two 12-core Xeon 2650-v4 CPUs, so a total of 24 cores (and as I noted, only 24 threads which avoids 2 separate threads fighting over access to a single FPU as there as only as many FPUs as there are cores).

I also ran it on a two 18core Xeon 2695v4 which has more cores and there is no performance gain. The performance is somewhat dominated by memory access (or cache misses) so you quickly hit the limit. When I run it on a two 26core Xeon Gold 6230R with 50% more memory bandwidth, I see the run-times for the parallel version drop by nearly 50%.

[Next page](https://chapel.discourse.group/t/code-optimization/17990.md?page=2)