/* Originally posted to povray.binaries.utilites by Grassblade, 13 Apr 2007 Version 0.02 The license of this file is POVray-ware, that is use of this file and/or any part of it is free, as long it is used directly in POVray or one of its unofficial rebuilds. Updated 2009.09.07 by CShake * fixed some places where #declare was used instead of #local * added colappend to the header comment, created rowappend * created Swap rows, Multiply row, Convert to Reduced Row Echelon Form, Gauss-Jordan Matrix Inversion Contents/ syntax: -Matricial sum syntax: summat(array1,array2,array3), effect: array3=array1+array2 -Scalar product syntax: scalprod(array1,array2,scalar) effect: array2=scalar*array1 -Transpose syntax: transp(array1,array2) effect: array2=array1' ('=transposed) -Matricial product syntax: matprod(array1, array2, array3) effect: array3=array1*array2 -Triple matricial product (useful for matrix inversion algorithm, henceforth UFMIL) syntax: triplematprod(array1,array2,array3,array4) effect: array4=((array1*array2)*array3) -Quadratic form syntax: quadform(vec1,array1,total) effect: total=(vec1)'*array1*vec1 -Row extract syntax: rowextr(array1,array2,rown) effect: array2 (row vector)=row number "rown" of array1, rown=1 means first row, etc. -Column extract syntax: colextr(array1,array2,coln) effect: array2 (column vector)=column number "coln" of array1, coln=1 means first column, etc. -Identity matrix syntax: identmat(matrixname, rank) effect: generates I(rank), a matrix filled with 1's on the main diagonal and 0's elsewhere -Matrix fill syntax: matrixfill(matrixname,m,n,value) effect: fills a m*n matrix ("matrixname") with constant "value" -Random fill: syntax: randomfill(matrixname,m,n,seedvalue) effect: fills a m*n matrix with random values generated with seed=seedvalue -Matrix subextract: syntax: subextract(array1,array2,a,b,c,d) effect: array2 is the matrix formed from extracting the values starting at row a, column b, to row c column d. a may be greater than c, and b of d. -Find maximum value (beta, may change): syntax: matrixmax(array1) effect: stores the highest value of the matrix in variable "maxi", also prints in message window the coordinates of the max -Find minimum value (beta, may change): syntax: matrixmin(array1) effect: stores the lowest value of the matrix in variable "mini", also prints in message window the coordinates of the min -ZZLY Special Bordered matrix (UFMIL) syntax: zzlyborder(array1,array2) effect: array2 is array1 plus an extra column and a row filled with 0's, except the intersection of the extra row and column has a value of 1 -Bordered Matrix: syntax: bordermatrix(array1,array2,array3,array4,array5) effect: array1 | array2 array5= -------+------- array3 | array4 ; #rows of array1 must equal #rows of array2, etc. -Matricial equality (UFMIL): syntax: matrixcond(array1,array2) effect: tests array1=array2 equality, the scalar "matrixcond1" will equal 1 if they are equal, 0 otherwise -Matrix (Column) append: syntax: colappend(array1,array2,array3) effect: array3= array1|array2 ; #rows must be the same for array1 and array2 -Matrix (Row) append: syntax: rowappend(array1,array2,array3) effect: array1 array3= ------ array2 ; #cols must be the same for array1 and array2 -Swap rows: syntax: swaprows(array1,r1,r2) effect: rows r1 and r2 are swapped in array1 -Multiply row: syntax: multiplyrow(array1,r,m) effect: multiplies every value in row r of array1 by m -Convert to Reduced Row Echelon Form: syntax: reducedrowechelon(array1,array2) effect: converts array1 to reduced row echelon form and stores result to array2 -Moore-Penrose generalised inverse: syntax: MPGinverse(array1,array2) effect: array2=(array1)^(-1) ; can invert singular and even rectangular matrices -Gauss-Jordan matrix inversion: syntax: GJinverse(array1,array2) effect: array2=(array1)^(-1) ; only works on square, nonsingular matrices note: considerably (~7.5x) faster than Moore-Penrose for matrices it can handle */ // matrix sum #macro summat(array1, array2,array3) #local dime=0; #local a1=dimension_size(array1,1); #if (a1=dimension_size(array2,1)) #local a2=dimension_size(array1,2); #if (a2=dimension_size(array2,2)) #local i=0; #declare array3= array[a1][a2]; #while (i<=a1-1) #local j=0; #while (j<=a2-1) #declare array3[i][j]=array1[i][j]+array2[i][j]; #local j=j+1; #end #local i=i+1; #end #else #local dime=1; #end #else #local dime=1; #end #if (dime=1) #debug "********Dimension mismatch (Sum)********\n" #end #end // scalar product #macro scalprod(array1,array2,scalar) #local a1=dimension_size(array1,1); #local a2=dimension_size(array1,2); #local i=0; #declare array2= array[a1][a2]; #while (i<=a1-1) #local j=0; #while (j<=a2-1) #declare array2[i][j]=array1[i][j]*scalar; #local j=j+1; #end #local i=i+1; #end #end // transpose #macro transp(array1,array2) #local a1=dimension_size(array1,1); #local a2=dimension_size(array1,2); #local i=0; #declare array2= array[a2][a1]; #while (i<=a2-1) #local j=0; #while (j<=a1-1) #declare array2[i][j]=array1[j][i]; #local j=j+1; #end #local i=i+1; #end #end // matricial product #macro matprod(array1, array2, array3) #local matm=dimension_size(array1,2); #if (matm=dimension_size(array2,1)) #local mato=dimension_size(array1,1); #local matn=dimension_size(array2,2); #local i=0; #declare array3= array[mato][matn]; #while (i<=mato-1) #local j=0; #while (j<=matn-1) #local summ=0; #local k=0; #while (k<=matm-1) #local summ=summ+array1[i][k]*array2[k][j]; #local k=k+1; #end #declare array3[i][j]=summ; #local j=j+1; #end #local i=i+1; #end #else #debug "********Dimension mismatch (Product)********\n" #end #end // triple product #macro triplematprod(array1,array2,array3,array4) #local matm=dimension_size(array1,2); #if (matm=dimension_size(array2,1)) matprod(array1,array2,temp1) #local matn=dimension_size(temp1,2); #if (matn=dimension_size(array3,1)) matprod(temp1,array3,array4) #else #debug "********Dimension mismatch (Triple Product)********\n" #end #else #debug "********Dimension mismatch (Triple Product)********\n" #end #end // quadratic forms vec1 is a column vector, array1 the square matrix, and total the name of the (scalar) result #macro quadform(vec1,array1,total) #local vecm=dimension_size(vec1,2); #if (vecm=1) #local vecn=dimension_size(vec1,1); #if (vecn=dimension_size(array1,1)) #local i=0; #declare total=0; #while (i<=vecn-1) #local j=0; #while (j<=vecn-1) #declare total=total+array1[i][j]*vec1[i][0]*vec1[j][0]; #local j=j+1; #end #local i=i+1; #end #else #debug "********Dimension mismatch (Quadratic Form)********\n" #end #else #debug "********First argument needs to be a vector**********\n" #end #end // row extract #macro rowextr(array1,array2,rown) // rown is an integer: 1=first row, etc. #local i=0; #declare array2=array[1][dimension_size(array1,2)]; #while (i<= dimension_size(array1,2)-1) #declare array2[0][i]=array1[rown-1][i]; #local i=i+1; #end #end /* column extract, array2 is the name of the column to be extracted from array1, coln is the ordinal of the column to be extracted: 1=first column, etc.*/ #macro colextr(array1,array2,coln) #local i=0; #declare array2=array[dimension_size(array1,1)][1]; #while (i<= dimension_size(array1,1)-1) #declare array2[i][0]=array1[i][coln-1]; #local i=i+1; #end #end // Identity matrix #macro identmat(matrixname, rank) #declare matrixname=array[rank][rank]; #local i=0; #while (i<= rank-1) #local j=0; #while (j<=rank-1) #if (i=j) #declare matrixname[i][j]=1; #else #declare matrixname[i][j]=0; #end #local j=j+1; #end #local i=i+1; #end #end // Fill matrix m rows, n columns with "value" #macro matrixfill(matrixname,m,n,value) #declare matrixname=array[m][n]; #local i=0; #while (i<= m-1) #local j=0; #while (j<=n-1) #declare matrixname[i][j]=value; #local j=j+1; #end #local i=i+1; #end #end // Random fill #macro randomfill(matrixname,m,n,seedvalue) #declare matrixname=array[m][n]; #declare rr5=seed(seedvalue); #local i=0; #while (i<= m-1) #local j=0; #while (j<=n-1) #declare matrixname[i][j]=rand(rr5); #local j=j+1; #end #local i=i+1; #end #end /* matrix subextract,array1 is starting matrix, array2 the name of the extracted matrix, a,b row and column start ordinals and c,d row and column end ordinals. If a > c they are swapped around, likewise for b and d.*/ #macro subextract(array1,array2,a,b,c,d) #if (a>c) #local e=a; #local a=c; #local c=e; #end #if (b>d) #local e=b; #local b=d; #local d=e; #end #declare array2=array[c-a+1][d-b+1]; #local i=0; #while (i<= c-a) #local j=0; #while (j<=d-b) #declare array2[i][j]=array1[a+i-1][b+j-1]; #local j=j+1; #end #local i=i+1; #end #end // max value #macro matrixmax(array1) #declare maxi=array1[0][0]; #local m=dimension_size(array1,1); #local n=dimension_size(array1,2); #local i=0; #while (i<= m-1) #local j=0; #while (j<=n-1) #if (array1[i][j]>maxi) #declare maxi=array1[i][j]; #local maxcoord1=i+1; #local maxcoord2=j+1; #end #local j=j+1; #end #local i=i+1; #end #debug concat("*********Matrix maximum: "str(maxi,1,0)" located at row "str(maxcoord1,-1,0)" and column "str(maxcoord2,-1,0)"*********\n" ) #end // min value #macro matrixmin(array1) #declare mini=array1[0][0]; #local m=dimension_size(array1,1); #local n=dimension_size(array1,2); #local i=0; #while (i<= m-1) #local j=0; #while (j<=n-1) #if (array1[i][j]array2 #macro matrixcond(array1,array2) #if (dimension_size(array1,1)=dimension_size(array2,1) & dimension_size(array1,2)=dimension_size(array2,2)) #local i=0; #declare matrixcond1=1; #local k=0; #while (i=0 & r2>=0 & r1=0 & r= numcol) #local row=numrow;// #debug "exiting1!\n" #else #local i = row; #local done=0; #while (atemp[i][lead] = 0 & done=0) #local i = i + 1; #if (i>=numrow) #local i = row; #local lead = lead + 1; #if (lead>=numcol) #local done=1; #local r=numrow;// #debug "exiting2!\n" #end //if #end //if #end //while #if(done=0) #if(i!=row) swaprows(atemp,i,row) #end #local mtemp=1/atemp[row][lead]; multiplyrow(atemp,row,mtemp) #local i=0; #while (i