matrix.h

fdtd 3D xyzPML MPI OpenMP

#include <stdlib.h> 
#include <malloc.h>

#ifndef __MatrixbyZsolt__

	#define __MatrixbyZsolt__

	///////////////////////////////////////////////////////////////////////////////////////
	//Initialize a 2D matrix
	///////////////////////////////////////////////////////////////////////////////////////
	template<class Tipus> Tipus ** Init_Matrix_2D(long  nx, long  ny)
	{
		Tipus **X;
		long  i;

		X = (Tipus **)calloc(nx,sizeof(Tipus *));
		if (X == NULL)
		{
			return NULL;
		}

		long  s = nx*ny; 

		Tipus* p = (Tipus*)calloc(s, sizeof(Tipus));
		if (p == NULL) 
		{
			free(X);
			return NULL;
		}
		
		for (i = 0; i < nx; i++)
		{
			X[i] = p;
			p += ny;
		}
		
		return X;
	};

	///////////////////////////////////////////////////////////////////////////////////////
	//Initialize a 2D matrix with different sizes in y
	///////////////////////////////////////////////////////////////////////////////////////
	template<class Tipus> Tipus ** Init_Matrix_2D(long  nx, long  *Size_y)
	{
	//the dimensions of Size_y must be [nx] 

		Tipus **X;
		long  i;

		X = (Tipus **)calloc(nx,sizeof(Tipus *));
		if (X == NULL)
		{
			return NULL;
		}

		long  s = 0;; 
		for (i = 0; i < nx; i++)
		{
			s += Size_y[i];
		}

		Tipus* p = (Tipus*)calloc(s, sizeof(Tipus));
		if (p == NULL) 
		{
			free(X);
			return NULL;
		}
		
		for (i = 0; i < nx; i++)
		{
			X[i] = p;
			p += Size_y[i];
		}
		
		return X;
	};

	///////////////////////////////////////////////////////////////////////////////////////
	//Clear a 2D matrix 
	///////////////////////////////////////////////////////////////////////////////////////
	template<class Tipus> Tipus ** Free_Matrix_2D(Tipus **X)
	{
		if (!X)
			return NULL;

		free(X[0]);
		free(X);
		return NULL;
	};
		
	///////////////////////////////////////////////////////////////////////////////////////
	//Initialize a 3D matrix
	///////////////////////////////////////////////////////////////////////////////////////
	template<class Tipus> Tipus *** Init_Matrix_3D(long  nx, long  ny, long  nz)
	{
		Tipus ***X;
		long  i, j, ii;

		X = (Tipus ***)calloc(nx,sizeof(Tipus **));
		if (X == NULL)
		{
			return NULL;
		}

		for (i = 0; i < nx; i++)
		{
			X[i] = (Tipus **)calloc(ny,sizeof(Tipus *));
			if (X[i] == NULL)
			{
				for (ii = 0; ii < i; ii++)
				{
					free(X[ii]);
				}
				free(X);
				return NULL;
			}
		}

		long  s = nx*ny*nz;
				
		Tipus* p = (Tipus*)calloc(s, sizeof(Tipus));
		if (p == NULL) 
		{
			for (i = 0; i< nx; i++)
			{
				free(X[i]);
			}
			free(X);
			return NULL;
		}
		
		for (i = 0; i < nx; i++)
		{
			for (j = 0; j < ny; j++)
			{
				X[i][j] = p;
				p += nz;
			}
		}

		return X;
	};

	///////////////////////////////////////////////////////////////////////////////////////
	//Initialize a 3D matrix with different sizes in z
	///////////////////////////////////////////////////////////////////////////////////////
	template<class Tipus> Tipus *** Init_Matrix_3D(long  nx, long  ny, long  ** Size_z)
	{
	//the dimensions of Size_z must be [nx][ny] 
		Tipus ***X;
		
		X = (Tipus ***)calloc(nx,sizeof(Tipus **));
		if (X==NULL)
		{
			return NULL;
		}

		long  i, j, ii;

		for (i = 0; i < nx; i++)
		{
			X[i] = (Tipus **)calloc(ny,sizeof(Tipus *));
			if (X[i] == NULL)
			{
				for (ii = 0; ii < i; ii++)
				{ 
					free(X[ii]);
				}
				free(X);
				return NULL;
			}
		}

		long  s = 0;
		for (i = 0; i < nx; i++)
		{
			for (j = 0; j < ny; j++)
			{
				s += Size_z[i][j];
			}
		}

		Tipus* p = (Tipus*)calloc(s, sizeof(Tipus));
		if (p == NULL) 
		{
			for (i = 0; i < nx; i++)
			{
				free(X[i]);
			}
			free(X);
			return NULL;
		}
		
		for (i = 0; i < nx; i++)
		{
			for (j = 0; j < ny; j++)
			{
				X[i][j] = p;
				p += Size_z[i][j];
			}
		}

		return X;
	};

	///////////////////////////////////////////////////////////////////////////////////////
	//Clear a 3D matrix 
	///////////////////////////////////////////////////////////////////////////////////////
	template<class Tipus> Tipus *** Free_Matrix_3D(Tipus ***X, long  nx)
	{
		if (!X)
			return NULL;


		free(X[0][0]);
		
		long  i;

		for (i = 0; i < nx; i++)
		{
			free(X[i]);
		}
 		free(X);

		return NULL;
	};

	///////////////////////////////////////////////////////////////////////////////////////
	//Initialize a 4D matrix
	///////////////////////////////////////////////////////////////////////////////////////
	template<class Tipus> Tipus **** Init_Matrix_4D(long  nx, long  ny, long  nz, long nw)
	{
		Tipus ****X;
		long  i, j, k, ii, jj, n_jj;

		//allocate in i
		X = (Tipus ****)calloc(nx,sizeof(Tipus ***));
		if (X == NULL)
		{
			return NULL;
		}

		//allocate in j
		for (i = 0; i < nx; i++)
		{
			X[i] = (Tipus ***)calloc(ny,sizeof(Tipus **));
			if (X[i] == NULL)
			{
				for (ii = 0; ii < i; ii++)
				{
					free(X[ii]);
				}
				free(X);
				return NULL;
			}
		}

		//allocate in k
		for (i = 0; i < nx; i++)
		{
			for (j = 0; j < ny; j++)
			{
				X[i][j] = (Tipus **)calloc(nz,sizeof(Tipus *));
				if (X[i][j] == NULL)
				{
					for (ii = 0; ii < i; ii++)
					{
						if (ii < i-1)
						{
							n_jj = ny;
						}
						else
						{
							n_jj = j;
						}

						for (jj = 0; jj < n_jj; jj++)
						{
							free(X[ii][jj]);
						}
						free(X[ii]);
					}
					free(X);
					return NULL;
				}
			}
		}

		long  s = nx*ny*nz*nw;
				
		Tipus* p = (Tipus*)calloc(s, sizeof(Tipus));
		if (p == NULL) 
		{
			for (i = 0; i < nx; i++)
			{
				for (j = 0; j < ny; j++)
				{
					free(X[i][j]);
				}
				free(X[i]);
			}
			free(X);
			return NULL;
		}
		
		for (i = 0; i < nx; i++)
		{
			for (j = 0; j < ny; j++)
			{
				for (k = 0; k < nz; k++)
				{
					X[i][j][k] = p;
					p += nw;
				}
			}
		}

		return X;
	};

	///////////////////////////////////////////////////////////////////////////////////////
	//Initialize a 4D matrix with different sizes in w
	///////////////////////////////////////////////////////////////////////////////////////
	template<class Tipus> Tipus **** Init_Matrix_4D(long  nx, long  ny, long  nz, long  ***Size_w)
	{
		//the dimensions of Size_w must be [nx][ny][nz] 
		Tipus ****X;
		
		X = (Tipus ****)calloc(nx,sizeof(Tipus ***));
		if (X == NULL)
		{
			return NULL;
		}

		long  i, j, k, ii, jj, n_jj;

		for (i = 0; i < nx; i++)
		{
			X[i] = (Tipus ***)calloc(ny,sizeof(Tipus **));
			if (X[i] == NULL)
			{
				for (ii = 0; ii < i; ii++)
				{ 
					free(X[ii]);
				}
				free(X);
				return NULL;
			}
		}

		for (i = 0; i < nx; i++)
		{
			for (j = 0; j < ny; j++)
			{
				X[i][j] = (Tipus **)calloc(nz,sizeof(Tipus *));
				if (X[i][j] == NULL)
				{
					for (ii = 0; ii < i; ii++)
					{
						if (ii < i-1)
						{
							n_jj = ny;
						}
						else
						{
							n_jj = j;
						}

						for (jj = 0; jj < n_jj; jj++)
						{
							free(X[ii][jj]);
						}
						free(X[ii]);
					}
					free(X);
					return NULL;
				}
			}
		}


		long  s = 0;
		for (i = 0; i < nx; i++)
		{
			for (j = 0; j < ny; j++)
			{
				for (k = 0; k < nz; k++)
				{
					s += Size_w[i][j][k];
				}
			}
		}

		Tipus* p = (Tipus*)calloc(s, sizeof(Tipus));
		if (p == NULL) 
		{
			for (i = 0; i < nx; i++)
			{
				for (j = 0; j < ny; j++)
				{
					free(X[i][j]);
				}
				free(X[i]);
			}
			free(X);
			return NULL;
		}
		
		for (i = 0; i < nx; i++)
		{
			for (j = 0; j < ny; j++)
			{
				for (k = 0; k < nz; k++)
				{
					X[i][j][k] = p;
					p += Size_w[i][j][k];
				}
			}
		}
		
		return X;
	};

	///////////////////////////////////////////////////////////////////////////////////////
	//Clear a 4D matrix 
	///////////////////////////////////////////////////////////////////////////////////////
	template<class Tipus> Tipus **** Free_Matrix_4D(Tipus ****X, long  nx, long  ny)
	{
		if (!X)
			return NULL;


		free(X[0][0][0]);
		
		long  i, j;

		for (i = 0; i < nx; i++)
		{
			for (j = 0; j < ny; j++)
			{
				free(X[i][j]);
			}
			free(X[i]);
		}
 		free(X);

		return NULL;
	};

#endif