 			_mx = mx[index+ number_of_contacts * 3];
 			_my = my[index+ number_of_contacts * 3];
 
-			_obj1 += _obj1_tmp * _mx;
+			obj1_temp[index+ number_of_contacts * 3] = _obj1_tmp * _mx;
 			_ms = _mx - _my;
 			ms[index + number_of_contacts * 3] = _ms;
 			_mg = mg[index+ number_of_contacts * 3];
-			temp1 += _mg * _ms;
-			temp2+= _ms * _ms;
+			temp1[index+ number_of_contacts * 3]= _mg * _ms;
+			temp2[index+ number_of_contacts * 3]= _ms * _ms;
 		}
 		if ( lm < min_val ) {
 #pragma critical

	}
	shurA(mx);
	real lm;
	real _obj1 = 0;
	real temp1 = 0;
	real temp2 = 0;

#pragma omp parallel private(lm)
	{
		lm=FLT_MAX;
#pragma omp for reduction(+:_obj1,temp1,temp2)
		for (uint index = 0; index < number_of_contacts; index++) {
			real3 temp = R3(0);
			int2 id_ = ids[index];


			real3 _obj1_tmp = 0.5 * _mg_tmp-_b;

			obj1_temp[index+ number_of_contacts * 0]= _obj1_tmp.x * _mx.x;
			obj1_temp[index+ number_of_contacts * 1]= _obj1_tmp.y * _mx.y;
			obj1_temp[index+ number_of_contacts * 2]= _obj1_tmp.z * _mx.z;

		}

#pragma omp for
		for (uint index = 0; index < number_of_bilaterals; index++) {
			real temp =0;
			int2 id_=ids[index+ number_of_contacts * 3];
