--- ray/src/common/mat4.c	1989/02/02 10:34:35	1.1
+++ ray/src/common/mat4.c	1990/01/18 23:58:56	1.4
@@ -13,7 +13,12 @@ static char SCCSid[] = "$SunId$ LBL";
 
 static double  m4tmp[4][4];		/* for efficiency */
 
+#ifdef  BSD
 #define  copymat4(m4a,m4b)	bcopy((char *)m4b,(char *)m4a,sizeof(m4tmp))
+#else
+#define  copymat4(m4a,m4b)	(void)memcpy((char *)m4a,(char *)m4b,sizeof(m4tmp))
+extern char  *memcpy();
+#endif
 
 
 setident4(m4)
@@ -51,13 +56,10 @@ double  v3a[3];
 register double  v3b[3];
 register double  m4[4][4];
 {
-	register int  i;
+	m4tmp[0][0] = v3b[0]*m4[0][0] + v3b[1]*m4[1][0] + v3b[2]*m4[2][0];
+	m4tmp[0][1] = v3b[0]*m4[0][1] + v3b[1]*m4[1][1] + v3b[2]*m4[2][1];
+	m4tmp[0][2] = v3b[0]*m4[0][2] + v3b[1]*m4[1][2] + v3b[2]*m4[2][2];
 	
-	for (i = 0; i < 3; i++)
-		m4tmp[0][i] = v3b[0]*m4[0][i] +
-			      v3b[1]*m4[1][i] +
-			      v3b[2]*m4[2][i];
-	
 	v3a[0] = m4tmp[0][0];
 	v3a[1] = m4tmp[0][1];
 	v3a[2] = m4tmp[0][2];
@@ -79,35 +81,37 @@ register double  m4[4][4];
 #ifdef  INVMAT
 /*
  * invmat - computes the inverse of mat into inverse.  Returns 1
- * if there exists an inverse, 0 otherwise.  It uses Gause Elimination
- * method.
+ * if there exists an inverse, 0 otherwise.  It uses Gaussian Elimination
+ * method with partial pivoting.
  */
 
 invmat(inverse,mat)
 double mat[4][4],inverse[4][4];
 {
 #define SWAP(a,b,t) (t=a,a=b,b=t)
+#define ABS(x) (x>=0?x:-(x))
 
 	register int i,j,k;
 	register double temp;
 
-	setident4(inverse);
 	copymat4(m4tmp, mat);
+	setident(inverse);
 
 	for(i = 0; i < 4; i++) {
-		if(m4tmp[i][i] == 0) {    /* Pivot is zero */
-			/* Look for a raw with pivot != 0 and swap raws */
-			for(j = i + 1; j < 4; j++)
-				if(m4tmp[j][i] != 0) {
-					for( k = 0; k < 4; k++) {
-						SWAP(m4tmp[i][k],m4tmp[j][k],temp);
-						SWAP(inverse[i][k],inverse[j][k],temp);
-						}
-					break;
-					}
-			if(j == 4)	/* No replacing raw -> no inverse */
-				return(0);
-			}
+		/* Look for row with largest pivot and swap rows */
+		temp = 0; j = -1;
+		for(k = i; k < 4; k++)
+			if(ABS(m4tmp[k][i]) > temp) {
+				temp = ABS(m4tmp[k][i]);
+				j = k;
+				}
+		if(j == -1)	/* No replacing row -> no inverse */
+			return(0);
+		if (j != i)
+			for(k = 0; k < 4; k++) {
+				SWAP(m4tmp[i][k],m4tmp[j][k],temp);
+				SWAP(inverse[i][k],inverse[j][k],temp);
+				}
 
 		temp = m4tmp[i][i];
 		for(k = 0; k < 4; k++) {
@@ -125,5 +129,8 @@ double mat[4][4],inverse[4][4];
 			}
 		}
 	return(1);
+
+#undef ABS
+#undef SWAP
 }
 #endif