""

svn path=/trunk/; revision=475

1 files changed, 202 insertions, 259 deletions

diff --git a/externals/grill/vasp/source/opdefs.h b/externals/grill/vasp/source/opdefs.h
index 18863f61..47006db6 100644
--- a/externals/grill/vasp/source/opdefs.h
+++ b/externals/grill/vasp/source/opdefs.h
@@ -20,345 +20,288 @@ WARRANTIES, see the file, "license.txt," in this distribution.
 #define _D_ALWAYS1 0
 #endif
 
+
 namespace VecOp {
 
 /*! \brief skeleton for unary real operations
 */
-template<class T,V FUN(T &v,T a)> BL _F__run(OpParam &p)										
+//template<class T,V FUN(T &v,T a)> BL V__run(register const T *sr,I rss,register T *dr,I rds,I frames)
+template<class T,class OP> BL V__run(register const T *sr,I rss,register T *dr,I rds,I frames)
 {																
-	register const S *sr = p.rsdt;								
-	register S *dr = p.rddt;									
-	register I i;												
-	if(sr == dr)												
-		if(_D_ALWAYS1 || p.rds == 1)							
-			_D_LOOP(i,p.frames)									
-            { FUN(*dr,*dr); dr++; }								
-			_E_LOOP												
+    register I i;
+    if(sr == dr)												
+		if(_D_ALWAYS1 || rds == 1)							
+            _DE_LOOP(i,frames, ( OP::run(*dr,*dr), dr++ ) )
 		else													
-			_D_LOOP(i,p.frames)									
-			{ FUN(*dr,*dr); dr += p.rds; }						
-			_E_LOOP												
+			_DE_LOOP(i,frames, ( OP::run(*dr,*dr), dr += rds ) )
 	else														
-		if(_D_ALWAYS1 || (p.rss == 1 && p.rds == 1))			
-			_D_LOOP(i,p.frames)									
-			{ FUN(*dr,*sr); sr++,dr++; }						
-			_E_LOOP												
+		if(_D_ALWAYS1 || (rss == 1 && rds == 1))			
+			_DE_LOOP(i,frames, ( OP::run(*dr,*sr), sr++,dr++ ) )
 		else													
-			_D_LOOP(i,p.frames)									
-			{ FUN(*dr,*sr); sr += p.rss,dr += p.rds; }			
-			_E_LOOP												
+			_DE_LOOP(i,frames, ( OP::run(*dr,*sr), sr += rss,dr += rds ) )
 	return true;												
 }
 
-template<class T,class CL> inline BL _D__run(OpParam &p) { return _F__run<T,CL::run>(p); }
-template<class T,class CL> inline BL d__run(OpParam &p) { return _d__run<T>(CL::run,p); }
 
 /*! \brief skeleton for unary complex operations
 */
-template<class T,V FUN(T &rv,T &iv,T ra,T ia)> BL _F__cun(OpParam &p)
+template<class T,class OP> BL V__cun(register const T *sr,register const T *si,I rss,I iss,register T *dr,register T *di,I rds,I ids,I frames)
 {																
-	register const S *sr = p.rsdt,*si = p.isdt;					
-	register S *dr = p.rddt,*di = p.iddt;						
-	register I i;												
+    register I i;
 	if(sr == dr && si == di)									
-		if(_D_ALWAYS1 || (p.rds == 1 && p.ids == 1))			
-			_D_LOOP(i,p.frames)									
-			{ FUN(*dr,*di,*dr,*di); dr++,di++; }				
-			_E_LOOP												
+		if(_D_ALWAYS1 || (rds == 1 && ids == 1))			
+            _DE_LOOP(i,frames, ( OP::cun(*dr,*di,*dr,*di), dr++,di++ ) )
 		else													
-			_D_LOOP(i,p.frames)									
-			{ FUN(*dr,*di,*dr,*di); dr += p.rds,di += p.ids; }	
-			_E_LOOP												
+			_DE_LOOP(i,frames, ( OP::cun(*dr,*di,*dr,*di), dr += rds,di += ids ) )
 	else														
-		if(_D_ALWAYS1 || (p.rss == 1 && p.iss == 1 && p.rds == 1 && p.ids == 1)) 
-			_D_LOOP(i,p.frames)									
-			{ FUN(*dr,*di,*sr,*si); sr++,si++,dr++,di++; }	
-			_E_LOOP												
+		if(_D_ALWAYS1 || (rss == 1 && iss == 1 && rds == 1 && ids == 1)) 
+			_DE_LOOP(i,frames, ( OP::cun(*dr,*di,*sr,*si), sr++,si++,dr++,di++ ) )
 		else													
-			_D_LOOP(i,p.frames)									
-			{ FUN(*dr,*di,*sr,*si);	sr += p.rss,si += p.iss,dr += p.rds,di += p.ids; } 
-			_E_LOOP												
+			_DE_LOOP(i,frames, ( OP::cun(*dr,*di,*sr,*si), sr += rss,si += iss,dr += rds,di += ids ) )
 	return true;												
 }
 
-template<class T,class CL> inline BL _D__cun(OpParam &p) { return _F__cun<T,CL::cun>(p); }
-template<class T,class CL> inline BL d__cun(OpParam &p) { return _d__cun<T>(CL::cun,p); }
+template<class T,class OP> BL V__vun(I layers,register const T *sr,register T *dr,I frames)
+{																
+    register I i;
+    switch(layers) {
+    case 1: 
+            V__run<T,OP>(sr,1,dr,1,frames); 
+            break;
+    case 2: 
+			_DF_LOOP(i,frames, ( vec_un<T,OP,2>(dr,sr), sr += 2, dr += 2) )
+            break;
+    case 3: 
+			_DF_LOOP(i,frames, ( vec_un<T,OP,3>(dr,sr), sr += 3, dr += 3) )
+            break;
+    case 4: 
+			_DF_LOOP(i,frames, ( vec_un<T,OP,4>(dr,sr), sr += 4, dr += 4) )
+            break;
+    default:
+			_DF_LOOP(i,frames, ( vec_un<T,OP,0>(dr,sr,layers), sr += layers, dr += layers) )
+            break;
+    }
+    return true;
+}
+
+
+template<class T,class OP> BL V__rbin(register const T *sr,I rss,register T *dr,I rds,register const T *ar,I ras,I frames)
+{																
+	register I i;												
+	if(sr == dr)								
+		if(_D_ALWAYS1 || (rds == 1 && ras == 1))				
+            _DE_LOOP(i,frames, ( OP::rbin(*dr,*dr,*ar), dr++,ar++ ) )
+		else												
+			_DE_LOOP(i,frames, ( OP::rbin(*dr,*dr,*ar), dr += rds,ar += ras ) )
+	else													
+		if(_D_ALWAYS1 || (rss == 1 && rds == 1 && ras == 1))	
+			_DE_LOOP(i,frames, ( OP::rbin(*dr,*sr,*ar), sr++,dr++,ar++ ) )
+		else												
+			_DE_LOOP(i,frames, ( OP::rbin(*dr,*sr,*ar), sr += rss,dr += rds,ar += ras ) )
+    return true;
+}
+
+template<class T,class OP> BL V__cbin(register const T *sr,register const T *si,I rss,I iss,register T *dr,register T *di,I ids,I rds,const T *ar,const T *ai,I ras,I ias,I frames)
+{																
+	register I i;												
+	if(sr == dr && si == di)							
+		if(_D_ALWAYS1 || (rds == 1 && ids == 1 && ras == 1 && ias == 1)) 
+            _DE_LOOP(i,frames, ( OP::cbin(*dr,*di,*dr,*di,*ar,*ai), dr++,di++,ar++,ai++ ) )
+		else											
+			_DE_LOOP(i,frames, ( OP::cbin(*dr,*di,*dr,*di,*ar,*ai), dr += rds,di += ids,ar += ras,ai += ias ) )
+	else												
+		_DE_LOOP(i,frames, ( OP::cbin(*dr,*di,*sr,*si,*ar,*ai), sr += rss,si += iss,dr += rds,di += ids,ar += ras,ai += ias ) )
+    return true;
+}
+
+
+template<class TR> class _A__vector {
+public:
+    static TR ev(const TR *a,I i,I m) { return a[i*m]; }
+};
+
+template<class TR> class _A__scalar {
+public:
+    static TR ev(TR a,I i,I m) { return a; }
+};
+
+class _A__env {
+public:
+    static R ev(Env::Iter &a,I i,I m) { return a.ValFwd(i); }
+};
+
+template<class T,class TA,class TR,class OP,class EVARG> BL Vx__rbin(register const T *sr,I rss,register T *dr,I rds,TA ar,I frames)
+{																
+	register I i;												
+	if(sr == dr)								
+		if(_D_ALWAYS1 || rds == 1)				
+            _DE_LOOP(i,frames, ( OP::rbin(*dr,*dr,EVARG::ev(ar,i,1)), dr++ ) )
+		else												
+            _DE_LOOP(i,frames, ( OP::rbin(*dr,*dr,EVARG::ev(ar,i,1)), dr += rds ) )
+	else													
+		if(_D_ALWAYS1 || rss == 1 && rds == 1)	
+            _DE_LOOP(i,frames, ( OP::rbin(*dr,*sr,EVARG::ev(ar,i,1)), sr++,dr++ ) )
+		else												
+            _DE_LOOP(i,frames, ( OP::rbin(*dr,*sr,EVARG::ev(ar,i,1)), sr += rss,dr += rds ) )
+    return true;
+}
+
+template<class T,class TA1,class TA2,class TR,class OP,class EVARG1,class EVARG2> BL Vx__cbin(register const T *sr,register const T *si,I rss,I iss,register T *dr,register T *di,I ids,I rds,TA1 ar,TA2 ai,I ras,I ias,I frames)
+{																
+	register I i;												
+	if(sr == dr && si == di)							
+		if(_D_ALWAYS1 || (rds == 1 && ids == 1 && ras == 1 && ias == 1)) 
+            _DE_LOOP(i,frames, ( OP::cbin(*dr,*di,*dr,*di,EVARG1::ev(ar,i,1),EVARG2::ev(ai,i,1)), dr++,di++ ) )
+		else											
+            _DE_LOOP(i,frames, ( OP::cbin(*dr,*di,*dr,*di,EVARG1::ev(ar,i,ras),EVARG2::ev(ai,i,ias)), dr += rds,di += ids ) )
+	else												
+        _DE_LOOP(i,frames, ( OP::cbin(*dr,*di,*sr,*si,EVARG1::ev(ar,i,ras),EVARG2::ev(ai,i,ias)), sr += rss,si += iss,dr += rds,di += ids ) )
+    return true;
+}
 
+template<class T,class TA,class TR,class OP,class EVARG> BL Vx__vbin(I layers,register const T *sr,register T *dr,TA ar,I frames)
+{																
+    register I i;
+    switch(layers) {
+    case 1: 
+            Vx__rbin<T,TA,TR,OP,EVARG>(sr,1,dr,1,ar,frames); 
+            break;
+    case 2: 
+			_DF_LOOP(i,frames, ( vec_bin<T,OP,2>(dr,sr,EVARG::ev(ar,i,2)), sr += 2, dr += 2) )
+            break;
+    case 3: 
+			_DF_LOOP(i,frames, ( vec_bin<T,OP,3>(dr,sr,EVARG::ev(ar,i,3)), sr += 3, dr += 3) )
+            break;
+    case 4: 
+			_DF_LOOP(i,frames, ( vec_bin<T,OP,4>(dr,sr,EVARG::ev(ar,i,4)), sr += 4, dr += 4) )
+            break;
+    default:
+			_DF_LOOP(i,frames, ( vec_bin<T,OP,0>(dr,sr,EVARG::ev(ar,i,layers),layers), sr += layers, dr += layers) )
+            break;
+    }
+    return true;
+}
+
+template<class T,class OP> inline BL V__vbin(I layers,register const T *sr,register T *dr,register const T *ar,I frames)
+{
+    return Vx__vbin<T,const T *,T,OP,_A__vector<T> >(layers,sr,dr,ar,frames);
+}
 
 /*! \brief skeleton for binary real operations
 */
-template<class T,V FUN(T &v,T a,T b)> BL _F__rbin(OpParam &p)					
+template<class T,class OP> BL _F__rbin(OpParam &p)					
 {																
-	register const S *sr = p.rsdt;								
-	register S *dr = p.rddt;									
-	register I i;												
 	if(p.HasArg() && p.arg[0].Is()) {							
 		switch(p.arg[0].argtp) {								
-		case OpParam::Arg::arg_v: {								
-			register const S *ar = p.arg[0].v.rdt;				
-			if(p.rsdt == p.rddt)								
-				if(_D_ALWAYS1 || (p.rds == 1 && p.arg[0].v.rs == 1))				
-					_D_LOOP(i,p.frames)									
-                    { FUN(*dr,*dr,*ar); dr++,ar++; }				
-					_E_LOOP												
-				else												
-					_D_LOOP(i,p.frames)									
-					{ FUN(*dr,*dr,*ar);	dr += p.rds,ar += p.arg[0].v.rs; } 
-					_E_LOOP												
-			else													
-				if(_D_ALWAYS1 || (p.rss == 1 && p.rds == 1 && p.arg[0].v.rs == 1))	
-					_D_LOOP(i,p.frames)									
-					{ FUN(*dr,*sr,*ar); sr++,dr++,ar++; }			
-					_E_LOOP												
-				else												
-					_D_LOOP(i,p.frames)									
-					{ FUN(*dr,*sr,*ar);	sr += p.rss,dr += p.rds,ar += p.arg[0].v.rs; } 
-					_E_LOOP												
+		case OpParam::Arg::arg_v: {	
+            V__rbin<T,OP>(p.rsdt,p.rss,p.rddt,p.rds,p.arg[0].v.rdt,p.arg[0].v.rs,p.frames);
 			break;												
 		}														
 		case OpParam::Arg::arg_env: {							
-			Env::Iter it(*p.arg[0].e.env); it.Init(0);			
-			if(p.rsdt == p.rddt)								
-				if(_D_ALWAYS1 || p.rds == 1)										
-					_D_LOOP(i,p.frames)									
-					{ FUN(*dr,*dr,it.ValFwd(i)); dr++; }			
-					_E_LOOP												
-				else												
-					_D_LOOP(i,p.frames)									
-					{ FUN(*dr,*dr,it.ValFwd(i)); dr += p.rds; }		
-					_E_LOOP												
-			else													
-				if(_D_ALWAYS1 || (p.rss == 1 && p.rds == 1))						
-					_D_LOOP(i,p.frames)									
-					{ FUN(*dr,*sr,it.ValFwd(i)); sr++,dr++; }		
-					_E_LOOP												
-				else												
-					_D_LOOP(i,p.frames)								 	
-					{ FUN(*dr,*sr,it.ValFwd(i)); sr += p.rss,dr += p.rds; }	
-					_E_LOOP												
+			Env::Iter it(*p.arg[0].e.env); it.Init(0);
+            Vx__rbin<T,Env::Iter &,R,OP,_A__env >(p.rsdt,p.rss,p.rddt,p.rds,it,p.frames);
 			break;												
 		}														
 		case OpParam::Arg::arg_x: {							
-			const R v =  p.arg[0].x.r;						
-			if(p.rsdt == p.rddt)							
-				if(_D_ALWAYS1 || p.rds == 1)				
-					_D_LOOP(i,p.frames)						
-					{ FUN(*dr,*dr,v); dr++; }			
-					_E_LOOP									
-				else										
-					_D_LOOP(i,p.frames)						
-					{ FUN(*dr,*dr,v); dr += p.rds; }	
-					_E_LOOP									
-			else											
-				if(_D_ALWAYS1 || (p.rss == 1 && p.rds == 1))
-					_D_LOOP(i,p.frames)						
-					{ FUN(*dr,*sr,v); sr++,dr++; }		
-					_E_LOOP									
-				else										
-					_D_LOOP(i,p.frames)						
-                    { FUN(*dr,*sr,v); sr += p.rss,dr += p.rds; }
-					_E_LOOP											
+            Vx__rbin<T,R,T,OP,_A__scalar<R> >(p.rsdt,p.rss,p.rddt,p.rds,p.arg[0].x.r,p.frames);
 			break;												
 		}														
 		}														
 	}															
 	else {														
-		register const S v = p.rbin.arg;						
-		if(p.rsdt == p.rddt)									
-			if(_D_ALWAYS1 || p.rds == 1)						
-				_D_LOOP(i,p.frames)									
-				{ FUN(*dr,*dr,v); dr++; }						
-				_E_LOOP												
-			else												
-				_D_LOOP(i,p.frames)									
-				{ FUN(*dr,*dr,v); dr += p.rds; }				
-				_E_LOOP												
-		else													
-			if(_D_ALWAYS1 || (p.rss == 1 && p.rds == 1))		
-				_D_LOOP(i,p.frames)									
-				{ FUN(*dr,*sr,v); sr++,dr++; }					
-				_E_LOOP												
-			else												
-				_D_LOOP(i,p.frames)									
-				{ FUN(*dr,*sr,v); sr += p.rss,dr += p.rds; }	
-				_E_LOOP												
+        Vx__rbin<T,T,T,OP,_A__scalar<T> >(p.rsdt,p.rss,p.rddt,p.rds,p.rbin.arg,p.frames);
 	}															
 	return true;												
 }
 
-template<class T,class CL> inline BL _D__rbin(OpParam &p) { return _F__rbin<T,CL::rbin>(p); }
-template<class T,class CL> inline BL d__rbin(OpParam &p) { return _d__rbin<T>(CL::rbin,p); }
-
 
 /*! \brief skeleton for binary complex operations
 */
-template<class T,V FUN(T &rv,T &iv,T ra,T ia,T rb,T ib)> BL _F__cbin(OpParam &p)							
+template<class T,class OP> BL _F__cbin(OpParam &p)							
 {																
-	register const S *sr = p.rsdt,*si = p.isdt;					
-	register S *dr = p.rddt,*di = p.iddt;						
-	register I i;												
 	if(p.HasArg() && p.arg[0].Is()) {							
 		switch(p.arg[0].argtp) {									
 		case OpParam::Arg::arg_v: {									
-			register const S *ar = p.arg[0].v.rdt,*ai = p.arg[0].v.idt;				
-			if(ai)													
-				if(sr == dr && si == di)							
-					if(_D_ALWAYS1 || (p.rds == 1 && p.ids == 1 && p.arg[0].v.rs == 1 && p.arg[0].v.is == 1)) 
-						_D_LOOP(i,p.frames)				 
-						{ FUN(*dr,*di,*dr,*di,*ar,*ai);	dr++,di++,ar++,ai++; }		
-						_E_LOOP												
-					else											
-						_D_LOOP(i,p.frames)				 
-						{ FUN(*dr,*di,*dr,*di,*ar,*ai);	dr += p.rds,di += p.ids,ar += p.arg[0].v.rs,ai += p.arg[0].v.is; }		
-						_E_LOOP												
-				else												
-					_D_LOOP(i,p.frames)				 
-					{ FUN(*dr,*di,*sr,*si,*ar,*ai);	sr += p.rss,si += p.iss,dr += p.rds,di += p.ids,ar += p.arg[0].v.rs,ai += p.arg[0].v.is; }			
-					_E_LOOP												
-			else													
-				if(sr == dr && si == di)							
-					_D_LOOP(i,p.frames)				 
-					{ FUN(*dr,*di,*dr,*di,*ar,0); dr += p.rds,di += p.ids,ar += p.arg[0].v.rs; }	
-					_E_LOOP												
-				else												
-					_D_LOOP(i,p.frames)				 
-					{ FUN(*dr,*di,*sr,*si,*ar,0); sr += p.rss,si += p.iss,dr += p.rds,di += p.ids,ar += p.arg[0].v.rs; } 
-					_E_LOOP												
+			if(p.arg[0].v.idt)
+                V__cbin<T,OP>(p.rsdt,p.isdt,p.rss,p.iss,p.rddt,p.iddt,p.rds,p.ids,p.arg[0].v.rdt,p.arg[0].v.idt,p.arg[0].v.rs,p.arg[0].v.is,p.frames);
+            else
+                Vx__cbin<T,const T *,T,T,OP,_A__vector<T>,_A__scalar<T> >(p.rsdt,p.isdt,p.rss,p.iss,p.rddt,p.iddt,p.rds,p.ids,p.arg[0].v.rdt,0,p.arg[0].v.rs,1,p.frames);
 			break;												
 		}														
 		case OpParam::Arg::arg_env: {							
 			Env::Iter it(*p.arg[0].e.env); it.Init(0);			
-			if(sr == dr && si == di)							
-				if(_D_ALWAYS1 || (p.rds == 1 && p.ids == 1)) 
-					_D_LOOP(i,p.frames)				 
-					{ FUN(*dr,*di,*dr,*di,it.ValFwd(i),0); dr++,di++; }			
-					_E_LOOP												
-				else											
-					_D_LOOP(i,p.frames)				 
-					{ FUN(*dr,*di,*dr,*di,it.ValFwd(i),0); dr += p.rds,di += p.ids; }			
-					_E_LOOP												
-			else												
-				_D_LOOP(i,p.frames)				 
-					{ FUN(*dr,*di,*sr,*si,it.ValFwd(i),0); sr += p.rss,si += p.iss,dr += p.rds,di += p.ids; }	
-				_E_LOOP												
+            Vx__cbin<T,Env::Iter &,T,R,OP,_A__env,_A__scalar<T> >(p.rsdt,p.isdt,p.rss,p.iss,p.rddt,p.iddt,p.rds,p.ids,it,0,1,1,p.frames);
 			break;												
 		}														
 		case OpParam::Arg::arg_x: {								
-			register const R ar = p.arg[0].x.r,ai = p.arg[0].x.i;				
-			if(sr == dr && si == di)							
-				if(_D_ALWAYS1 || (p.rds == 1 && p.ids == 1)) 
-					_D_LOOP(i,p.frames)				 
-					{ FUN(*dr,*di,*dr,*di,ar,ai); dr++,di++; }			
-					_E_LOOP												
-				else											
-					_D_LOOP(i,p.frames)				 
-					{ FUN(*dr,*di,*dr,*di,ar,ai); dr += p.rds,di += p.ids; }			
-					_E_LOOP												
-			else												
-				_D_LOOP(i,p.frames)				 
-				{ FUN(*dr,*di,*sr,*si,ar,ai); sr += p.rss,si += p.iss,dr += p.rds,di += p.ids; }	
-				_E_LOOP												
+            Vx__cbin<T,R,R,R,OP,_A__scalar<R>,_A__scalar<R> >(p.rsdt,p.isdt,p.rss,p.iss,p.rddt,p.iddt,p.rds,p.ids,p.arg[0].x.r,p.arg[0].x.i,1,1,p.frames);
 			break;												
 		}														
 		}														
 	}															
 	else {														
-		register const S rv = p.cbin.rarg,iv = p.cbin.iarg;		
-		if(sr == dr && si == di)								
-			if(_D_ALWAYS1 || (p.rds == 1 && p.ids == 1))		
-				_D_LOOP(i,p.frames)						
-				{ FUN(*dr,*di,*dr,*di,rv,iv); dr++,di++; }					
-				_E_LOOP												
-			else												
-				_D_LOOP(i,p.frames)				 
-				{ FUN(*dr,*di,*dr,*di,rv,iv); dr += p.rds,di += p.ids; }					
-				_E_LOOP												
-		else													
-			if(_D_ALWAYS1 || (p.rds == 1 && p.ids == 1 && p.rss == 1 && p.iss == 1)) 
-				_D_LOOP(i,p.frames)				 
-				{ FUN(*dr,*di,*sr,*si,rv,iv); sr++,si++,dr++,di++; }					
-				_E_LOOP												
-			else												
-				_D_LOOP(i,p.frames)				 
-				{ FUN(*dr,*di,*sr,*si,rv,iv); sr += p.rss,si += p.iss,dr += p.rds,di += p.ids; }					
-				_E_LOOP												
+        Vx__cbin<T,T,T,T,OP,_A__scalar<T>,_A__scalar<T> >(p.rsdt,p.isdt,p.rss,p.iss,p.rddt,p.iddt,p.rds,p.ids,p.cbin.rarg,p.cbin.iarg,1,1,p.frames);
 	}															
 	return true;												
 }
 
-template<class T,class CL> inline BL _D__cbin(OpParam &p) { return _F__cbin<T,CL::cbin>(p); }			
-template<class T,class CL> inline BL d__cbin(OpParam &p) { return _d__cbin<T>(CL::cbin,p); }
-
-
 /*! \brief skeleton for real operations with parameter block
 */
-template<class T,V FUN(T &v,T r,OpParam &p)> BL _F__rop(OpParam &p)						
+template<class T,class ARG,class OP> BL V__rop(ARG p,register const S *sr,I rss,register S *dr,I rds,I frames)
 {																
-	register const S *sr = p.rsdt;								
-	register S *dr = p.rddt;									
 	register I i;												
 	if(sr == dr)												
-		if(_D_ALWAYS1 || p.rds == 1)											
-			_D_LOOP(i,p.frames)										
-            { FUN(*dr,*dr,p); dr++; }							
-			_E_LOOP												
+		if(_D_ALWAYS1 || rds == 1)											
+            _DE_LOOP(i,frames, ( OP::rop(*dr,*dr,p), dr++ ) )
 		else													
-			_D_LOOP(i,p.frames)									
-			{ FUN(*dr,*dr,p); dr += p.rds; }					
-			_E_LOOP												
+			_DE_LOOP(i,frames, ( OP::rop(*dr,*dr,p), dr += rds ) )
 	else														
-		if(_D_ALWAYS1 || (p.rss == 1 && p.rds == 1))			
-			_D_LOOP(i,p.frames)									
-			{ FUN(*dr,*sr,p); sr++,dr++; }					
-			_E_LOOP												
+		if(_D_ALWAYS1 || (rss == 1 && p.rds == 1))			
+			_DE_LOOP(i,frames, ( OP::rop(*dr,*sr,p), sr++,dr++ ) )
 		else													
-			_D_LOOP(i,p.frames)									
-			{ FUN(*dr,*sr,p); sr += p.rss,dr += p.rds; }		
-			_E_LOOP												
+			_DE_LOOP(i,frames, ( OP::rop(*dr,*sr,p), sr += rss,dr += rds ) )
 	return true;												
 }
 
-template<class T,class CL> inline BL _D__rop(OpParam &p) { return _F__rop<T,CL::rop>(p); }			
-template<class T,class CL> inline BL d__rop(OpParam &p) { return _d__rop<T>(CL::rop,p); }
-
-
 /*! \brief skeleton for complex operations with parameter block
 */
-template<class T,V FUN(T &rv,T &iv,T ra,T ia,OpParam &p)> BL _F__cop(OpParam &p)						
+template<class T,class ARG,class OP> BL V__cop(ARG p,register const S *sr,register const S *si,I rss,I iss,register S *dr,register S *di,I rds,I ids,I frames)
 {																
-	register const S *sr = p.rsdt,*si = p.isdt;					
-	register S *dr = p.rddt,*di = p.iddt;						
 	register I i;												
 	if(sr == dr && si == di)									
-		if(_D_ALWAYS1 || (p.rds == 1 && p.ids == 1))			
-			_D_LOOP(i,p.frames)									
-			{ FUN(*dr,*di,*dr,*di,p); dr++,di++; }			
-			_E_LOOP												
+		if(_D_ALWAYS1 || (rds == 1 && ids == 1))			
+			_DE_LOOP(i,frames, ( OP::cop(*dr,*di,*dr,*di,p), dr++,di++ ) )
 		else													
-			_D_LOOP(i,p.frames)									
-			{ FUN(*dr,*di,*dr,*di,p); dr += p.rds,di += p.ids; } 
-			_E_LOOP												
+			_DE_LOOP(i,frames, ( OP::cop(*dr,*di,*dr,*di,p), dr += rds,di += ids ) )
 	else														
 		if(_D_ALWAYS1 || (p.rss == 1 && p.iss == 1 && p.rds == 1 && p.ids == 1)) 
-			_D_LOOP(i,p.frames)				 
-			{ FUN(*dr,*di,*sr,*si,p); sr++,si++,dr++,di++; }	
-			_E_LOOP												
+			_DE_LOOP(i,frames, ( OP::cop(*dr,*di,*sr,*si,p), sr++,si++,dr++,di++ ) )
 		else													
-			_D_LOOP(i,p.frames)				 
-			{ FUN(*dr,*di,*sr,*si,p); sr += p.rss,si += p.iss,dr += p.rds,di += p.ids; }	
-			_E_LOOP												
+			_DE_LOOP(i,frames, ( OP::cop(*dr,*di,*sr,*si,p), sr += rss,si += iss,dr += rds,di += ids ) )
 	return true;												
 }
 
-template<class T,class CL> inline BL _D__cop(OpParam &p) { return _F__cop<T,CL::cop>(p); }
+template<class T,class CL> inline BL _D__run(OpParam &p) { return V__run<T,CL>(p.rsdt,p.rss,p.rddt,p.rds,p.frames); }
+template<class T,class CL> inline BL _D__cun(OpParam &p) { return V__cun<T,CL>(p.rsdt,p.isdt,p.rss,p.iss,p.rddt,p.iddt,p.rds,p.ids,p.frames); }
+template<class T,class CL> inline BL _D__rbin(OpParam &p) { return _F__rbin<T,CL>(p); }
+template<class T,class CL> inline BL _D__cbin(OpParam &p) { return _F__cbin<T,CL>(p); }			
+template<class T,class CL> inline BL _D__rop(OpParam &p) { return V__rop<T,OpParam &,CL>(p,p.rsdt,p.rss,p.rddt,p.rds,p.frames); }
+template<class T,class CL> inline BL _D__cop(OpParam &p) { return V__cop<T,OpParam &,CL>(p,p.rsdt,p.isdt,p.rss,p.iss,p.rddt,p.iddt,p.rds,p.ids,p.frames); }
+
+
+template<class T,class CL> inline BL d__run(OpParam &p) { return _d__run<T>(CL::run,p); }
+template<class T,class CL> inline BL d__cun(OpParam &p) { return _d__cun<T>(CL::cun,p); }
+template<class T,class CL> inline BL d__rbin(OpParam &p) { return _d__rbin<T>(CL::rbin,p); }
+template<class T,class CL> inline BL d__cbin(OpParam &p) { return _d__cbin<T>(CL::cbin,p); }
+template<class T,class CL> inline BL d__rop(OpParam &p) { return _d__rop<T>(CL::rop,p); }
 template<class T,class CL> inline BL d__cop(OpParam &p) { return _d__cop<T>(CL::cop,p); }
 
-template<class T,V FUN(T &v,T a)> inline BL f__run(OpParam &p) { return _d__run(FUN,p); }
-template<class T,V FUN(T &rv,T &iv,T ra,T ia)> inline BL f__cun(OpParam &p) { return _d__cun(FUN,p); }
-template<class T,V FUN(T &v,T a,T b)> inline BL f__rbin(OpParam &p) { return _d__rbin(FUN,p); }
-template<class T,V FUN(T &rv,T &iv,T ra,T ia,T rb,T ib)> inline BL f__cbin(OpParam &p) { return _d__cbin(FUN,p); }
-template<class T,V FUN(T &v,T r,OpParam &p)> inline BL f__rop(OpParam &p) { return _d__rop(FUN,p); }
-template<class T,V FUN(T &rv,T &iv,T ra,T ia,OpParam &p)> inline BL f__cop(OpParam &p) { return _d__cop(FUN,p); }
+template<class T,class CL> inline BL f__run(OpParam &p) { return _d__run(CL::run,p); }
+template<class T,class CL> inline BL f__cun(OpParam &p) { return _d__cun(CL::cun,p); }
+template<class T,class CL> inline BL f__rbin(OpParam &p) { return _d__rbin(CL::rbin,p); }
+template<class T,class CL> inline BL f__cbin(OpParam &p) { return _d__cbin(CL::cbin,p); }
+template<class T,class CL> inline BL f__rop(OpParam &p) { return _d__rop(CL::rop,p); }
+template<class T,class CL> inline BL f__cop(OpParam &p) { return _d__cop(CL::cop,p); }
 
 
 #ifdef VASP_COMPACT
@@ -368,12 +311,12 @@ template<class T,class CL> inline BL D__rbin(OpParam &p) { return _d__rbin<T>(CL
 template<class T,class CL> inline BL D__cbin(OpParam &p) { return _d__cbin<T>(CL::cbin,p); }
 template<class T,class CL> inline BL D__rop(OpParam &p) { return _d__rop<T>(CL::rop,p); }
 template<class T,class CL> inline BL D__cop(OpParam &p) { return _d__cop<T>(CL::cop,p); }
-template<class T,V FUN(T &v,T a)> inline BL F__run(OpParam &p) { return _d__run(FUN,p); }
-template<class T,V FUN(T &rv,T &iv,T ra,T ia)> inline BL F__cun(OpParam &p) { return _d__cun<T>(FUN,p); }
-template<class T,V FUN(T &v,T a,T b)> inline BL F__rbin(OpParam &p) { return _d__rbin<T>(FUN,p); }
-template<class T,V FUN(T &rv,T &iv,T ra,T ia,T rb,T ib)> inline BL F__cbin(OpParam &p) { return _d__cbin<T>(FUN,p); }
-template<class T,V FUN(T &v,T r,OpParam &p)> inline BL F__rop(OpParam &p) { return _d__rop<T>(FUN,p); }
-template<class T,V FUN(T &rv,T &iv,T ra,T ia,OpParam &p)> inline BL F__cop(OpParam &p) { return _d__cop<T>(FUN,p); }
+template<class T,class CL> inline BL F__run(OpParam &p) { return _d__run(CL::run,p); }
+template<class T,class CL> inline BL F__cun(OpParam &p) { return _d__cun<T>(CL::cun,p); }
+template<class T,class CL> inline BL F__rbin(OpParam &p) { return _d__rbin<T>(CL::rbin,p); }
+template<class T,class CL> inline BL F__cbin(OpParam &p) { return _d__cbin<T>(CL::cbin,p); }
+template<class T,class CL> inline BL F__rop(OpParam &p) { return _d__rop<T>(CL::rop,p); }
+template<class T,class CL> inline BL F__cop(OpParam &p) { return _d__cop<T>(CL::cop,p); }
 #else
 template<class T,class CL> inline BL D__run(OpParam &p) { return _D__run<T,CL>(p); }
 template<class T,class CL> inline BL D__cun(OpParam &p) { return _D__cun<T,CL>(p); }
@@ -381,12 +324,12 @@ template<class T,class CL> inline BL D__rbin(OpParam &p) { return _D__rbin<T,CL>
 template<class T,class CL> inline BL D__cbin(OpParam &p) { return _D__cbin<T,CL>(p); }
 template<class T,class CL> inline BL D__rop(OpParam &p) { return _D__rop<T,CL>(p); }
 template<class T,class CL> inline BL D__cop(OpParam &p) { return _D__cop<T,CL>(p); }
-template<class T,V FUN(T &v,T a)> inline BL F__run(OpParam &p) { return _F__run<T,FUN>(p); }
-template<class T,V FUN(T &rv,T &iv,T ra,T ia)> inline BL F__cun(OpParam &p) { return _F__cun<T,FUN>(p); }
-template<class T,V FUN(T &v,T a,T b)> inline BL F__rbin(OpParam &p) { return _F__rbin<T,FUN>(p); }
-template<class T,V FUN(T &rv,T &iv,T ra,T ia,T rb,T ib)> inline BL F__cbin(OpParam &p) { return _F__cbin<T,FUN>(p); }
-template<class T,V FUN(T &v,T r,OpParam &p)> inline BL F__rop(OpParam &p) { return _F__rop<T,FUN>(p); }
-template<class T,V FUN(T &rv,T &iv,T ra,T ia,OpParam &p)> inline BL F__cop(OpParam &p) { return _F__cop<T,FUN>(p); }
+template<class T,class CL> inline BL F__run(OpParam &p) { return V__run<T,CL>(p.rsdt,p.rss,p.rddt,p.rds,p.frames); }
+template<class T,class CL> inline BL F__cun(OpParam &p) { return V__cun<T,CL>(p.rsdt,p.isdt,p.rss,p.iss,p.rddt,p.iddt,p.rds,p.ids,p.frames); }
+template<class T,class CL> inline BL F__rbin(OpParam &p) { return _F__rbin<T,CL>(p); }
+template<class T,class CL> inline BL V__cbin(OpParam &p) { return _F__cbin<T,CL>(p); }
+template<class T,class CL> inline BL F__rop(OpParam &p) { return V__rop<T,OpParam &,CL>(p,p.rsdt,p.rss,p.rddt,p.rds,p.frames); }
+template<class T,class CL> inline BL F__cop(OpParam &p) { return V__cop<T,OpParam &,CL>(p,p.rsdt,p.isdt,p.rss,p.iss,p.rddt,p.iddt,p.rds,p.ids,p.frames); }
 #endif
 
 } // namespace VecOp