Xref: utzoo comp.sources.wanted:12734 sci.math:11883 Path: utzoo!utgpu!news-server.csri.toronto.edu!cs.utexas.edu!yale!cmcl2!lanl!nmsu!opus!ted From: ted@nmsu.edu (Ted Dunning) Newsgroups: comp.sources.wanted,sci.math Subject: Re: Inverse of 4x4 Matrix Message-ID: Date: 3 Aug 90 17:39:33 GMT References: <1990Aug3.021200.21023@oswego.Oswego.EDU> Sender: news@NMSU.Edu Followup-To: comp.sources.wanted Organization: NMSU Computer Science Lines: 180 In-reply-to: cs-ind04@penelope.Oswego.EDU's message of 3 Aug 90 02:12:00 GMT From: cs-ind04@penelope.Oswego.EDU (Jason E. Forbes) Newsgroups: comp.sources.wanted,sci.math Before I reinvent the wheel, does anyone have a C function that can find the inverse of a 4 by 4 matrix? ... the real answer to your question is to buy yourself a copy of numerical recipes in c. just for fun, though, i revved maple up and had it spit out the inversion in fortran and then used emacs to change it into c. the total time including cursory debugging was about 30 minutes, and the resulting code obviously would provide a major opportunity for a fancy optimizer, but would probably be a disaster on a vectorizer. but, of course, 4x4 inversion isn't what you would be doing on a vector machine anyway. no guarantees for numerical stability, but anything at all should work pretty well for such small matrices. it would be an interesting question to to an operation count on the resulting compiled code and comparing the result to what a looping version would do. it would also be interesting to compare the number of operations that a common subexpression optimizer could save. so anyway, have fun with this. /* invert a 4x4 matrix stored in a and store the result in b. both a and b must be preallocated. this code was generated automatically in fortran and then translated to c, so it should be thoroughly tested before use */ invert_4_4(a,b) double **a, **b; { double t1, t2, t4, t5, t6, t7, t8, t9, t11, t12, t14, t15, t17; double t19, t22, t23, t26, t27, t29, t30, t32, t33, t48, t50, t53, t56; double t58, t59, t61, t62, t66, t69, t71, t73, t76, t79, t84, t116, t118; double t121, t124, t126, t128, t133, t135, t138, t141, t143, t145, t151; double t153, t159, t161, t166, t168, t174, t176, t187, t189, t212, t214; t1 = a[2][2]*a[3][3]; t2 = a[0][1]*t1; t4 = a[2][3]*a[3][2]; t5 = a[0][1]*t4; t6 = a[0][2]*a[3][3]; t7 = a[2][1]*t6; t8 = a[0][3]*a[3][2]; t9 = a[2][1]*t8; t11 = a[0][2]*a[2][3]; t12 = a[3][1]*t11; t14 = a[0][3]*a[2][2]; t15 = a[3][1]*t14; t17 = a[1][1]*t1; t19 = a[1][1]*t4; t22 = a[1][2]*a[3][3]; t23 = a[2][1]*t22; t26 = a[1][3]*a[3][2]; t27 = a[2][1]*t26; t29 = a[1][2]*a[2][3]; t30 = a[3][1]*t29; t32 = a[1][3]*a[2][2]; t33 = a[3][1]*t32; t48 = a[0][1]*t22; t50 = a[0][1]*t26; t53 = a[1][1]*t6; t56 = a[1][1]*t8; t58 = a[0][2]*a[1][3]; t59 = a[3][1]*t58; t61 = a[0][3]*a[1][2]; t62 = a[3][1]*t61; t66 = a[0][1]*t29; t69 = a[0][1]*t32; t71 = a[1][1]*t11; t73 = a[1][1]*t14; t76 = a[2][1]*t58; t79 = a[2][1]*t61; t84 = 1/(a[0][0]*t17-a[0][0]*t19-a[0][0]*t23+a[0][0]*t27+a[0][0]*t30-a[0][0]*t33-a[1][0]*t2+a[1][0]*t5+a[1][0]*t7-a[1][0]*t9-a[1][0]*t12+a[1][0]*t15+a[2][0]*t48-a[2][0]*t50-a[2][0]*t53+a[2][0]*t56+a[2][0]*t59-a[2][0]*t62-a[3][0]*t66+a[3][0]*t69+a[3][0]*t71-a[3][0]*t73-a[3][0]*t76+a[3][0]*t79); t116 = a[2][1]*a[3][2]; t118 = a[2][2]*a[3][1]; t121 = a[0][1]*a[3][2]; t124 = a[0][2]*a[3][1]; t126 = a[0][1]*a[2][2]; t128 = a[0][2]*a[2][1]; t133 = a[2][1]*a[3][3]; t135 = a[2][3]*a[3][1]; t138 = a[0][1]*a[3][3]; t141 = a[0][3]*a[3][1]; t143 = a[0][1]*a[2][3]; t145 = a[0][3]*a[2][1]; t151 = a[1][1]*a[2][2]; t153 = a[1][2]*a[2][1]; t159 = a[0][1]*a[1][2]; t161 = a[0][2]*a[1][1]; t166 = a[1][1]*a[3][3]; t168 = a[1][3]*a[3][1]; t174 = a[0][1]*a[1][3]; t176 = a[0][3]*a[1][1]; t187 = a[1][1]*a[2][3]; t189 = a[1][3]*a[2][1]; t212 = a[1][1]*a[3][2]; t214 = a[1][2]*a[3][1]; b[0][1] = (-t2+t5+t7-t9-t12+t15)*t84; b[0][2] = -(-t48+t50+t53-t56-t59+t62)*t84; b[1][2] = -(a[0][0]*t22-a[0][0]*t26-a[1][0]*t6+a[1][0]*t8+a[3][0]*t58-a[3][0]*t61)*t84; b[1][3] = -(-a[0][0]*t29+a[0][0]*t32+a[1][0]*t11-a[1][0]*t14-a[2][0]*t58+a[2][0]*t61)*t84; b[3][1] = (a[0][0]*t116-a[0][0]*t118-a[2][0]*t121+a[2][0]*t124+a[3][0]*t126-a[3][0]*t128)*t84; b[2][1] = -(a[0][0]*t133-a[0][0]*t135-a[2][0]*t138+a[2][0]*t141+a[3][0]*t143-a[3][0]*t145)*t84; b[3][3] = (a[0][0]*t151-a[0][0]*t153-a[1][0]*t126+a[1][0]*t128+a[2][0]*t159-a[2][0]*t161)*t84; b[2][2] = (a[0][0]*t166-a[0][0]*t168-a[1][0]*t138+a[1][0]*t141+a[3][0]*t174-a[3][0]*t176)*t84; b[2][0] = (a[1][0]*t133-a[1][0]*t135-a[2][0]*t166+a[2][0]*t168+a[3][0]*t187-a[3][0]*t189)*t84; b[1][0] = -(a[1][0]*t1-a[1][0]*t4-a[2][0]*t22+a[2][0]*t26+a[3][0]*t29-a[3][0]*t32)*t84; b[0][3] = -(t66-t69-t71+t73+t76-t79)*t84; b[3][2] = -(a[0][0]*t212-a[0][0]*t214-a[1][0]*t121+a[1][0]*t124+a[3][0]*t159-a[3][0]*t161)*t84; b[0][0] = (t17-t19-t23+t27+t30-t33)*t84; b[2][3] = -(a[0][0]*t187-a[0][0]*t189-a[1][0]*t143+a[1][0]*t145+a[2][0]*t174-a[2][0]*t176)*t84; b[3][0] = -(a[1][0]*t116-a[1][0]*t118-a[2][0]*t212+a[2][0]*t214+a[3][0]*t151-a[3][0]*t153)*t84; b[1][1] = (a[0][0]*t1-a[0][0]*t4-a[2][0]*t6+a[2][0]*t8+a[3][0]*t11-a[3][0]*t14)*t84; } /* multiply a times b to get c. c must be preallocated. code generated in fortran by maple and then converted semi-automatically to c. it may well be buggy. */ mul_4_4(a,b,c) double **a, **b, **c; { c[0][1] = a[0][0]*b[0][1]+a[0][1]*b[1][1]+a[0][2]*b[2][1]+a[0][3]*b[3][1]; c[0][2] = a[0][0]*b[0][2]+a[0][1]*b[1][2]+a[0][2]*b[2][2]+a[0][3]*b[3][2]; c[1][2] = a[1][0]*b[0][2]+a[1][1]*b[1][2]+a[1][2]*b[2][2]+a[1][3]*b[3][2]; c[1][3] = a[1][0]*b[0][3]+a[1][1]*b[1][3]+a[1][2]*b[2][3]+a[1][3]*b[3][3]; c[3][1] = a[3][0]*b[0][1]+a[3][1]*b[1][1]+a[3][2]*b[2][1]+a[3][3]*b[3][1]; c[2][1] = a[2][0]*b[0][1]+a[2][1]*b[1][1]+a[2][2]*b[2][1]+a[2][3]*b[3][1]; c[3][3] = a[3][0]*b[0][3]+a[3][1]*b[1][3]+a[3][2]*b[2][3]+a[3][3]*b[3][3]; c[2][2] = a[2][0]*b[0][2]+a[2][1]*b[1][2]+a[2][2]*b[2][2]+a[2][3]*b[3][2]; c[2][0] = a[2][0]*b[0][0]+a[2][1]*b[1][0]+a[2][2]*b[2][0]+a[2][3]*b[3][0]; c[1][0] = a[1][0]*b[0][0]+a[1][1]*b[1][0]+a[1][2]*b[2][0]+a[1][3]*b[3][0]; c[0][3] = a[0][0]*b[0][3]+a[0][1]*b[1][3]+a[0][2]*b[2][3]+a[0][3]*b[3][3]; c[3][2] = a[3][0]*b[0][2]+a[3][1]*b[1][2]+a[3][2]*b[2][2]+a[3][3]*b[3][2]; c[0][0] = a[0][0]*b[0][0]+a[0][1]*b[1][0]+a[0][2]*b[2][0]+a[0][3]*b[3][0]; c[2][3] = a[2][0]*b[0][3]+a[2][1]*b[1][3]+a[2][2]*b[2][3]+a[2][3]*b[3][3]; c[3][0] = a[3][0]*b[0][0]+a[3][1]*b[1][0]+a[3][2]*b[2][0]+a[3][3]*b[3][0]; c[1][1] = a[1][0]*b[0][1]+a[1][1]*b[1][1]+a[1][2]*b[2][1]+a[1][3]*b[3][1]; } /* simple cursory test of the inversion and multiplication routines */ main() { double **a,**b,**c; void *calloc(); int i,j; a = calloc(4,sizeof(a[0])); b = calloc(4,sizeof(b[0])); c = calloc(4,sizeof(c[0])); for (i=0;i<4;i++) { a[i] = calloc(4,sizeof(a[0][0])); b[i] = calloc(4,sizeof(b[0][0])); c[i] = calloc(4,sizeof(c[0][0])); for (j=0;j<4;j++) { a[i][j] = random()%1000; } } invert_4_4(a,b); mul_4_4(a,b,c); for (i=0;i<4;i++) { for (j=0;j<4;j++) { printf("%19.17lf ", c[i][j]); } printf("\n"); } } -- Offer void except where prohibited by law.