/* Copyright Henrik Stewenius, 2005. All Rights Reserved. Permission to use [*] this software for academic purposes is hereby granted without fee or royalty, subject to the following conditions. 1. The software is provided "as is" by the copyright holder, with absolutely no actual or implied warranty of correctness, fitness, intellectual property ownership, or anything else whatsoever. You use the software entirely at your own risk. In no event shall the copyright holder be liable for any direct, indirect or perceived damage whatsoever connected with the use of this work. 2. The software may not be distributed [*], whether in its original or modified form, without explicit written consent from the copyright holder. 3. The copyright holder retains all copyright to the work. 4. Any scholarly work, research or publication making use of this software must explicitly acknowledge the copyright holder. In the case of a publication, citation of the work, and author thereof, titled: @INPROCEEDINGS{stewenius-etal-omnivis-2005, AUTHOR = {H. Stew\'enius and D. Nist\'er and M. Oskarsson and K. {\AA}str{\"o}m}, TITLE = {Solutions to Minimal Generalized Relative Pose Problems}, BOOKTITLE = {Workshop on Omnidirectional Vision}, YEAR = 2005, MONTH = OCT, ADDRESS = {Beijing China}, PDF = {http://www.vis.uky.edu/~stewe/publications/stewenius_05_omnivis_sm26gen.pdf} } and/or @PHDTHESIS{stewenius-pdh-2005, AUTHOR = {H. Stew\'enius}, TITLE = {Gr{\"o}bner Basis Methods for Minimal Problems in Computer Vision}, SCHOOL = {Lund University}, YEAR = 2005, MONTH = APR, PDF = {http://www.maths.lth.se/matematiklth/personal/stewe/THESIS/stewenius_phd_2005.pdf} } is deemed sufficient acknowledgement. Permission to use or distribute [*] this software for commercial or military purposes, or under any conditions not covered by the conditions 1-4 above, is NOT granted. Such permission may be obtained subject to explicit prior written consent from the copyright holder. [*] `Use' denotes any use whatsoever including downloading, compilation, execution and modification. `Distribution' includes any distribution of the software and any derived works such as object files, executable programs and libraries. */ #include "mex.h" #include #define K000 9 #define K001 8 #define K002 5 #define K010 7 #define K011 4 #define K020 2 #define K100 6 #define K101 3 #define K110 1 #define K200 0 #define KK21_000 KK[ n1_30+(1-1)*10+K000] #define KK21_001 KK[ n1_30+(1-1)*10+K001] #define KK21_002 KK[ n1_30+(1-1)*10+K002] #define KK21_010 KK[ n1_30+(1-1)*10+K010] #define KK21_011 KK[ n1_30+(1-1)*10+K011] #define KK21_020 KK[ n1_30+(1-1)*10+K020] #define KK21_100 KK[ n1_30+(1-1)*10+K100] #define KK21_101 KK[ n1_30+(1-1)*10+K101] #define KK21_110 KK[ n1_30+(1-1)*10+K110] #define KK21_200 KK[ n1_30+(1-1)*10+K200] #define KK22_000 KK[ n1_30+(2-1)*10+K000] #define KK22_001 KK[ n1_30+(2-1)*10+K001] #define KK22_002 KK[ n1_30+(2-1)*10+K002] #define KK22_010 KK[ n1_30+(2-1)*10+K010] #define KK22_011 KK[ n1_30+(2-1)*10+K011] #define KK22_020 KK[ n1_30+(2-1)*10+K020] #define KK22_100 KK[ n1_30+(2-1)*10+K100] #define KK22_101 KK[ n1_30+(2-1)*10+K101] #define KK22_110 KK[ n1_30+(2-1)*10+K110] #define KK22_200 KK[ n1_30+(2-1)*10+K200] #define KK23_000 KK[ n1_30+(3-1)*10+K000] #define KK23_001 KK[ n1_30+(3-1)*10+K001] #define KK23_002 KK[ n1_30+(3-1)*10+K002] #define KK23_010 KK[ n1_30+(3-1)*10+K010] #define KK23_011 KK[ n1_30+(3-1)*10+K011] #define KK23_020 KK[ n1_30+(3-1)*10+K020] #define KK23_100 KK[ n1_30+(3-1)*10+K100] #define KK23_101 KK[ n1_30+(3-1)*10+K101] #define KK23_110 KK[ n1_30+(3-1)*10+K110] #define KK23_200 KK[ n1_30+(3-1)*10+K200] #define KK31_000 KK[ n2_30+(1-1)*10+K000] #define KK31_001 KK[ n2_30+(1-1)*10+K001] #define KK31_002 KK[ n2_30+(1-1)*10+K002] #define KK31_010 KK[ n2_30+(1-1)*10+K010] #define KK31_011 KK[ n2_30+(1-1)*10+K011] #define KK31_020 KK[ n2_30+(1-1)*10+K020] #define KK31_100 KK[ n2_30+(1-1)*10+K100] #define KK31_101 KK[ n2_30+(1-1)*10+K101] #define KK31_110 KK[ n2_30+(1-1)*10+K110] #define KK31_200 KK[ n2_30+(1-1)*10+K200] #define KK32_000 KK[ n2_30+(2-1)*10+K000] #define KK32_001 KK[ n2_30+(2-1)*10+K001] #define KK32_002 KK[ n2_30+(2-1)*10+K002] #define KK32_010 KK[ n2_30+(2-1)*10+K010] #define KK32_011 KK[ n2_30+(2-1)*10+K011] #define KK32_020 KK[ n2_30+(2-1)*10+K020] #define KK32_100 KK[ n2_30+(2-1)*10+K100] #define KK32_101 KK[ n2_30+(2-1)*10+K101] #define KK32_110 KK[ n2_30+(2-1)*10+K110] #define KK32_200 KK[ n2_30+(2-1)*10+K200] #define KK33_000 KK[ n2_30+(3-1)*10+K000] #define KK33_001 KK[ n2_30+(3-1)*10+K001] #define KK33_002 KK[ n2_30+(3-1)*10+K002] #define KK33_010 KK[ n2_30+(3-1)*10+K010] #define KK33_011 KK[ n2_30+(3-1)*10+K011] #define KK33_020 KK[ n2_30+(3-1)*10+K020] #define KK33_100 KK[ n2_30+(3-1)*10+K100] #define KK33_101 KK[ n2_30+(3-1)*10+K101] #define KK33_110 KK[ n2_30+(3-1)*10+K110] #define KK33_200 KK[ n2_30+(3-1)*10+K200] #define KK41_000 KK[ n3_30+(1-1)*10+K000] #define KK41_001 KK[ n3_30+(1-1)*10+K001] #define KK41_002 KK[ n3_30+(1-1)*10+K002] #define KK41_010 KK[ n3_30+(1-1)*10+K010] #define KK41_011 KK[ n3_30+(1-1)*10+K011] #define KK41_020 KK[ n3_30+(1-1)*10+K020] #define KK41_100 KK[ n3_30+(1-1)*10+K100] #define KK41_101 KK[ n3_30+(1-1)*10+K101] #define KK41_110 KK[ n3_30+(1-1)*10+K110] #define KK41_200 KK[ n3_30+(1-1)*10+K200] #define KK42_000 KK[ n3_30+(2-1)*10+K000] #define KK42_001 KK[ n3_30+(2-1)*10+K001] #define KK42_002 KK[ n3_30+(2-1)*10+K002] #define KK42_010 KK[ n3_30+(2-1)*10+K010] #define KK42_011 KK[ n3_30+(2-1)*10+K011] #define KK42_020 KK[ n3_30+(2-1)*10+K020] #define KK42_100 KK[ n3_30+(2-1)*10+K100] #define KK42_101 KK[ n3_30+(2-1)*10+K101] #define KK42_110 KK[ n3_30+(2-1)*10+K110] #define KK42_200 KK[ n3_30+(2-1)*10+K200] #define KK43_000 KK[ n3_30+(3-1)*10+K000] #define KK43_001 KK[ n3_30+(3-1)*10+K001] #define KK43_002 KK[ n3_30+(3-1)*10+K002] #define KK43_010 KK[ n3_30+(3-1)*10+K010] #define KK43_011 KK[ n3_30+(3-1)*10+K011] #define KK43_020 KK[ n3_30+(3-1)*10+K020] #define KK43_100 KK[ n3_30+(3-1)*10+K100] #define KK43_101 KK[ n3_30+(3-1)*10+K101] #define KK43_110 KK[ n3_30+(3-1)*10+K110] #define KK43_200 KK[ n3_30+(3-1)*10+K200] int V1[]={0,1,2,3,4,5,6,8,9,10,11,12,13,15,16,17,18,19,21,22,23,24,26,27,28,30,31,33,35,36,37,38,39,40,42,43,44,45,46,48,49,50,51,53,54,55,58,59,61,64,65,66,67,68,70,71,72,73,75,76,77,79,80,82,85,86,87,88,90,91,92,94,95,97,100,101,102,104,105,107,110,111,113,116}; int V2[]={1,2,3,4,5,6,7,9,10,11,12,13,14,16,17,18,19,20,22,23,24,25,27,28,29,31,32,34,36,37,38,39,40,41,43,44,45,46,47,49,50,51,52,54,55,57,59,60,62,65,66,67,68,69,71,72,73,74,76,77,78,80,81,83,86,87,88,89,91,92,93,95,96,98,101,102,103,105,106,108,111,112,114,117}; int V3[]={8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,56,42,43,44,45,46,47,48,49,50,51,52,53,54,55,57,58,59,60,61,62,63,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,90,91,92,93,94,95,96,97,98,99,104,105,106,107,108,109,113,114,115,118}; int It[]={35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,57,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87,88,89,90,91,92,93,94,95,96,97,98,99,100,101,102,103,104,105,106,107,108,109,110,111,112,113,114,115,116,117,118,119}; int ETT[]={35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,57,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87,88,89,90,91,92,93,94,95,96,97,98,99,100,101,102,103,104,105,106,107,108,109,110,111,112,113,114,115,116,117,118,119}; int REORDER[]={0,1,2,3,4,5,6,8,9,10,11,12,13,15,16,22,23,24,25,26,27,7,28,29,30,31,32,14,33,17,34,35,36,37,38,39,40,41,42,43,44,18,19,20,21,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59}; int N1[] ={2-1, 2-1, 2-1, 3-1, 2-1, 2-1, 3-1, 2-1, 3-1, 4-1} ; int N2[] ={3-1, 3-1, 4-1, 4-1, 3-1, 4-1, 4-1, 5-1, 5-1, 5-1} ; int N3[] ={4-1, 5-1, 5-1, 5-1, 6-1, 6-1, 6-1, 6-1, 6-1, 6-1} ; int fixed_block[]={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }; /* The above block was built using the following matlab code: fprintf('int V1[]={'); for i=1:84, fprintf('%i,', (find(V1(i,I))) -1), end,fprintf('};\n'); fprintf('int V2[]={'); for i=1:84, fprintf('%i,', (find(V2(i,I))) -1), end,fprintf('};\n'); fprintf('int V3[]={'); for i=1:84, fprintf('%i,', (find(V3(i,I))) -1), end,fprintf('};\n'); fprintf('int ETT[]={'); for i=1:84, fprintf('%i,', (find(ETT(i,I))) -1), end,fprintf('};\n'); II1 = [1:7 22 8:13 28 14 15 30 42:45] ; II2 =setdiff( 1:60 , II1 ); II = [II1 II2]; [dummy, iii ] = sort( II); fprintf('int REORDER[]={'); fprintf('%i,', iii-1); fprintf('}\n\n') */ mxArray *elim1( mxArray *AA); void mult_U_w_all( double *Unew, double *B, int Um); void elim2( double *B); void buildA( double *A, double *L1, double *L2, int Nblocks ); void build_last_gb( double *restrow, double *U, int N); void build_first_row_of_M( double *restrow, double *U ,double *last_gb, int N); #ifdef MAIN_BUILD_EQ2 void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { double *L1, *L2, *A; mxArray *AA; int Nblocks=2; if(nrhs <1 ) mexErrMsgTxt("L\n"); if(nrhs >2 ) mexErrMsgTxt("L,Nblocks\n"); if(nrhs >1 ) Nblocks = mxGetScalar( prhs[1] ); L1 = mxGetPr(mxGetCell( prhs[0] , 0 )) ; L2 = mxGetPr(mxGetCell( prhs[0] , 1 )) ; AA = mxCreateDoubleMatrix(10*Nblocks ,84 , mxREAL); A = mxGetPr( AA); buildA( A,L1, L2 , Nblocks) ; plhs[0] =AA; } #endif #ifdef MAIN_MULT_WITH_ALL void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { double *U,*B; int N,m; if(nrhs !=1 ) mexErrMsgTxt("U\n"); m= mxGetM(prhs[0]); /* if( mxGetM(prhs[0]) != 15 ) mexErrMsgTxt("H%G�%@den skall vara 15.\n");*/ if( mxGetN(prhs[0]) != 84 ) mexErrMsgTxt("Bredden skall vara 84.\n"); U = mxGetPr(prhs[0]); plhs[0] = mxCreateDoubleMatrix( m*4,120, mxREAL); B = mxGetPr(plhs[0]); mult_U_w_all( U, B, m); } #endif #ifdef MAIN_INITM void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { double *M; plhs[0] = mxCreateDoubleMatrix( 64,64, mxREAL); M = mxGetPr( plhs[0] ) ; M[154]=1; M[219]=1; M[348]=1; M[541]=1; M[606]=1; M[671]=1; M[736]=1; M[801]=1; M[930]=1; M[995]=1; M[1060]=1; M[1125]=1; M[1254]=1; M[1319]=1; M[1384]=1; M[1513]=1; M[1578]=1; M[1707]=1; M[1900]=1; M[1965]=1; M[2030]=1; M[2095]=1; M[2224]=1; M[2289]=1; M[2354]=1; M[2483]=1; M[2548]=1; M[2677]=1; M[2870]=1; M[2935]=1; M[3000]=1; M[3129]=1; M[3194]=1; M[3323]=1; M[3516]=1; M[3581]=1; M[3710]=1; M[3903]=1; } #endif #ifdef MAIN_LASTGB void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { mxArray *UU; double *U, *last_gb; int N; plhs[0] = mxCreateDoubleMatrix( 1,120, mxREAL); last_gb = mxGetPr( plhs[0] ) ; if(nrhs !=1 ) mexErrMsgTxt("U\n"); UU = prhs[0]; U = mxGetPr( UU ); N = mxGetM( UU); build_last_gb(last_gb, U, N); } #endif #ifdef MAIN_FIRST_ROW void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { int N; mxArray *UU; double *U, *last_gb,*restrow; if(nrhs !=2 ) mexErrMsgTxt("U,last_gb\n"); UU = prhs[0]; U = mxGetPr( UU ); N = mxGetM( UU); last_gb = mxGetPr( prhs[1] ); plhs[0] = mxCreateDoubleMatrix( 1,120, mxREAL); restrow = mxGetPr( plhs[0]); build_first_row_of_M( restrow, U ,last_gb,N); } #endif mxArray *elim1( mxArray *AA){ double *A,*U, *Unew,*B; mxArray *invec[2],*L, *outvec[2], *UU, *UUnew, *BB; double t, temp, max_val, max_val_signed; int m,n,row, i,j, max_I; A = mxGetPr( AA ); m = mxGetM( AA ); n = mxGetN( AA ); if( n != 84 ) mexErrMsgTxt("A skall vara 84 bred!\n"); if( m%10 ) mexErrMsgTxt("A skall ha höjd delbar 10!\n"); if( m < 30 ) mexErrMsgTxt("A skall ha minst 30 rader!\n"); outvec[0] = AA; mexCallMATLAB( 2, invec, 1, outvec,"lu" ); L = invec[0]; UU = invec[1]; mxDestroyArray( L ); U = mxGetPr( UU); for(row=0; row < 15 ; row++ ){ /* Först sätta förstaelementet till 1*/ t = U[row + 30*row]; for(i=0; i< n; i++){ U[row + 30*i] = U[row + 30*i]/t; } /*Nu skall det göras radoperationer */ for(j=0; j< row; j++){ t = U[j + 30*row]; for(i=0; i< n; i++){ U[j + 30*i] -= t*U[row + 30*i]; } } } return UU; } void mult_U_w_all( double *Unew, double *B, int Um){ int i,j,Um2,Um3,Um4; int temp_reorderETT; int temp_reorderV1; int temp_reorderV2; int temp_reorderV3; double current; Um2 = 2*Um; Um3 = 3*Um; Um4 = 4*Um; /*Slänger in alla multiplikationerna i en slinga. Se koden nedan för en mera begriplig version.*/ for(i=0; i< Um; i++){ temp_reorderETT = REORDER[Um3+i]; temp_reorderV3 = REORDER[(Um2)+i]; temp_reorderV2 = REORDER[Um+i]; temp_reorderV1 = REORDER[i]; for( j=0; j< 84; j++){ current = Unew[i+Um*j]; B[temp_reorderETT+Um4*ETT[j]] = current; B[temp_reorderV1+Um4*V1[j]] = current; B[temp_reorderV2+Um4*V2[j]] = current; B[temp_reorderV3+Um4*V3[j]] = current; } } } void mult_U_w_all_old( double *Unew, double *B, int Um){ int i,j,Um2,Um3,Um4, temp_reorder; Um2 = 2*Um; Um3 = 3*Um; Um4 = 4*Um; /*Det borde gå att skriva lite snabbare genom att slå ihop de 4 multiplikationerna till en slinga istället. */ /*Tar den lätta multiplikation med ETT först */ for(i=0; i< Um; i++){ temp_reorder = REORDER[Um3+i]; for( j=0; j< 84; j++){ B[temp_reorder+(Um4)*ETT[j]] = Unew[i+Um*j]; } } /*Nu skall det multipliceras med V1 */ for(i=0; i< Um; i++){ temp_reorder = REORDER[i]; for( j=0; j< 84; j++){ B[temp_reorder+(Um4)*V1[j]] = Unew[i+Um*j]; } } /*Nu skall det multipliceras med V2 */ for(i=0; i< Um; i++){ temp_reorder = REORDER[Um+i]; for( j=0; j< 84; j++){ B[temp_reorder+(Um4)*V2[j]] = Unew[i+Um*j]; } } /*Nu skall det multipliceras med V3 */ for(i=0; i< Um; i++){ temp_reorder = REORDER[(Um2)+i]; for( j=0; j< 84; j++){ B[temp_reorder+(Um4)*V3[j]] = Unew[i+Um*j]; } } } void elim2( double *B){ int row, i, j , max_I; double t, max_val, max_val_signed, temp; /*Radoperationer */ for(row=1; row < 7; row++){ for(i=0; i< 120; i++){ B[21+row + 60*i] -= B[row + 60*i]; } } for(row=9; row < 14; row++){ for(i=0; i< 120; i++){ B[19+row + 60*i] -= B[row + 60*i]; } } for(row=8; row < 17; row++){ for(i=0; i< 120; i++){ B[26+row + 60*i] -= B[row + 60*i]; } } for(i=0; i< 120; i++){ B[33 + 60*i] -= B[16 + 60*i]; } /*Nu är vi klara med den triviala delen av eliminationen! */ /*Nu är det dags att inleda arbetet! */ for( row = 17; row < 56; row ++){ if( (row >= 35) && (row<50) ){ max_I = row+10; max_val=1; max_val_signed =1; }else{ max_val=-1; for(i=row; i< 60; i++){ if( fabs( B[i + row*60] ) > max_val ){ max_I = i; max_val = fabs( B[i + row*60]); } } max_val_signed = B[max_I + row*60]; } if( max_val<0) mexErrMsgTxt("pivotering sket sig!"); if( row != max_I){ for(i=0; i< 120; i++){ temp = B[row + 60*i]; B[row + 60*i]= B[max_I + 60*i]/max_val_signed; B[max_I + 60*i] = temp; } }else{ for(i=0; i< 120; i++){ B[row + 60*i]= B[row + 60*i]/max_val_signed; } } /*Pivotering klar. Dags att snurra. */ for(j=0; j< 60; j++){ if( row!= j ){ t = B[j + 60*row]; for(i=0; i< 120; i++){ B[j + 60*i] -= t* B[row + 60*i]; } } } } } #if 1 void buildA( double *A, double *L1, double *L2 , int Nblocks){ double KK[6*3*10]; int row,i1,i2,i3,n1_30,n2_30,n3_30, cnt, fixed_point,blockid; double a11, b11, a12, b12;/* a13, b13, a14, b14, a15, b15, a16, b16;*/ double a21, b21, a22, b22;/* a23, b23, a24, b24, a25, b25, a26, b26;*/ double a31, b31, a32, b32;/* a33, b33, a34, b34, a35, b35, a36, b36;*/ double a41, b41, a42, b42;/* a43, b43, a44, b44, a45, b45, a46, b46;*/ double a51, b51, a52, b52;/* a53, b53, a54, b54, a55, b55, a56, b56;*/ double a61, b61, a62, b62;/* a63, b63, a64, b64, a65, b65, a66, b66;*/ for(blockid = 0 ; blockid < Nblocks; blockid++){ fixed_point = blockid; row = 0; for( cnt = 0; cnt <= 5 ; cnt++ ){ if( cnt != fixed_point ){ row++; a11 = L1[0+6*fixed_point]; b11 = L2[0+6*fixed_point]; a12 = L1[0+6*cnt]; b12 = L2[0+6*cnt]; a21 = L1[1+6*fixed_point]; b21 = L2[1+6*fixed_point]; a22 = L1[1+6*cnt]; b22 = L2[1+6*cnt]; a31 = L1[2+6*fixed_point]; b31 = L2[2+6*fixed_point]; a32 = L1[2+6*cnt]; b32 = L2[2+6*cnt]; a41 = L1[3+6*fixed_point]; b41 = L2[3+6*fixed_point]; a42 = L1[3+6*cnt]; b42 = L2[3+6*cnt]; a51 = L1[4+6*fixed_point]; b51 = L2[4+6*fixed_point]; a52 = L1[4+6*cnt]; b52 = L2[4+6*cnt]; a61 = L1[5+6*fixed_point]; b61 = L2[5+6*fixed_point]; a62 = L1[5+6*cnt]; b62 = L2[5+6*cnt]; KK[row*30+0*10 + K200 ] = -b22*a52-b32*a62-b62*a32-b22*a32*b31*b51+b12*a32*b11*b61-b12*a32*b31*b41-b32*a22*b21*b61+b32*a12*b11*b61+b42*a12-b12*a22*a21*a41-b12*a32*a31*a41+b12*a32*a11*a61+b12*a22*a11*a51-b32*a12*b31*b41+b32*a22*b31*b51-b52*a22+b12*a42+b12*a22*b11*b51-b12*a22*b21*b41+b22*a32*b21*b61-b22*a12*b21*b41+b22*a12*b11*b51-b22*a12*a21*a41+b22*a12*a11*a51+b22*a32*a31*a51-b22*a32*a21*a61-b32*a22*a31*a51+b32*a22*a21*a61-b32*a12*a31*a41+b32*a12*a11*a61; KK[row*30+0*10 + K110 ] = 2*b12*a12*b21*b41-2*b12*a12*b11*b51-2*b12*a12*a11*a51-2*b12*a32*a31*a51+2*b12*a32*a21*a61-2*b32*a12*b31*b51+2*b12*a12*a21*a41-2*b22*a22*a21*a41+2*b22*a22*a11*a51-2*b22*a32*a31*a41+2*b32*a12*b21*b61-2*b32*a22*b31*b41+2*b32*a22*b11*b61+2*b22*a22*b11*b51-2*b22*a22*b21*b41+2*b22*a32*a11*a61+2*b52*a12+2*b42*a22+2*b22*a42+2*b12*a52; KK[row*30+0*10 + K020 ] = b22*a52-b32*a62-b62*a32-b22*a32*b31*b51+b12*a32*b11*b61-b12*a32*b31*b41+b32*a22*b21*b61-b32*a12*b11*b61-b42*a12+b12*a22*a21*a41+b12*a32*a31*a41-b12*a32*a11*a61-b12*a22*a11*a51+b32*a12*b31*b41-b32*a22*b31*b51+b52*a22-b12*a42-b12*a22*b11*b51+b12*a22*b21*b41+b22*a32*b21*b61+b22*a12*b21*b41-b22*a12*b11*b51+b22*a12*a21*a41-b22*a12*a11*a51-b22*a32*a31*a51+b22*a32*a21*a61-b32*a22*a31*a51+b32*a22*a21*a61-b32*a12*a31*a41+b32*a12*a11*a61; KK[row*30+0*10 + K101 ] = 2*b12*a12*a31*a41-2*b12*a12*a11*a61+2*b12*a22*a31*a51-2*b12*a22*a21*a61+2*b22*a12*b31*b51-2*b12*a12*b11*b61+2*b32*a22*a11*a51-2*b22*a32*b21*b41+2*b22*a32*b11*b51-2*b22*a12*b21*b61-2*b32*a32*b31*b41+2*b32*a32*b11*b61-2*b32*a32*a31*a41+2*b32*a32*a11*a61-2*b32*a22*a21*a41+2*b12*a12*b31*b41+2*b62*a12+2*b42*a32+2*b32*a42+2*b12*a62; KK[row*30+0*10 + K011 ] = 2*b32*a52+2*b22*a12*a31*a41+2*b22*a22*a31*a51-2*b22*a12*a11*a61-2*b12*a32*b11*b51-2*b32*a12*a11*a51+2*b32*a12*a21*a41-2*b22*a22*a21*a61+2*b32*a32*a21*a61-2*b32*a32*a31*a51-2*b12*a22*b11*b61+2*b12*a32*b21*b41+2*b32*a32*b21*b61+2*b12*a22*b31*b41-2*b22*a22*b21*b61+2*b22*a22*b31*b51-2*b32*a32*b31*b51+2*b62*a22+2*b22*a62+2*b52*a32; KK[row*30+0*10 + K002 ] = -b42*a12+b32*a22*b31*b51+b32*a12*b31*b41-b32*a12*b11*b61-b32*a22*b21*b61-b12*a22*a11*a51-b12*a32*a11*a61+b12*a32*a31*a41+b12*a22*a21*a41+b32*a62-b12*a22*b21*b41+b12*a22*b11*b51+b12*a32*b31*b41-b52*a22-b12*a32*b11*b61+b62*a32+b22*a32*b31*b51-b22*a12*b11*b51+b22*a12*b21*b41-b22*a32*b21*b61-b12*a42-b22*a52+b32*a12*a31*a41-b32*a22*a21*a61+b32*a22*a31*a51-b22*a32*a21*a61+b22*a32*a31*a51-b32*a12*a11*a61-b22*a12*a21*a41+b22*a12*a11*a51; KK[row*30+0*10 + K100 ] = -2*b12*a32*b11*b51+2*b12*a32*b21*b41+2*b12*a22*b11*b61-2*b22*a12*a11*a61-2*b12*a22*b31*b41-2*b22*a22*b31*b51-2*b32*a32*b31*b51+2*b32*a32*b21*b61-2*b32*a12*a21*a41+2*b32*a12*a11*a51+2*b32*a32*a31*a51-2*b32*a32*a21*a61+2*b22*a22*a31*a51-2*b22*a22*a21*a61+2*b22*a22*b21*b61+2*b22*a12*a31*a41+2*b52*a32+2*b22*a62-2*b62*a22-2*b32*a52; KK[row*30+0*10 + K010 ] = -2*b12*a62+2*b32*a42-2*b42*a32+2*b62*a12-2*b22*a12*b21*b61+2*b12*a12*a11*a61-2*b12*a12*a31*a41+2*b12*a12*b31*b41-2*b32*a22*a21*a41+2*b32*a32*a11*a61-2*b32*a32*a31*a41-2*b32*a32*b11*b61+2*b32*a32*b31*b41-2*b22*a32*b11*b51+2*b22*a32*b21*b41-2*b12*a12*b11*b61+2*b22*a12*b31*b51+2*b32*a22*a11*a51+2*b12*a22*a21*a61-2*b12*a22*a31*a51; KK[row*30+0*10 + K001 ] = -2*b22*a42-2*b52*a12+2*b42*a22+2*b12*a52+2*b12*a32*a21*a61-2*b12*a32*a31*a51-2*b12*a12*a11*a51+2*b12*a12*a21*a41-2*b32*a12*b21*b61+2*b32*a12*b31*b51+2*b22*a22*b11*b51-2*b22*a22*b21*b41+2*b32*a22*b11*b61-2*b32*a22*b31*b41-2*b22*a32*a11*a61+2*b22*a32*a31*a41-2*b22*a22*a11*a51+2*b22*a22*a21*a41+2*b12*a12*b11*b51-2*b12*a12*b21*b41; KK[row*30+0*10 + K000 ] = b42*a12-b32*a22*b31*b51-b32*a12*b31*b41+b32*a12*b11*b61+b32*a22*b21*b61+b12*a22*a11*a51+b12*a32*a11*a61-b12*a32*a31*a41-b12*a22*a21*a41+b32*a62+b12*a22*b21*b41-b12*a22*b11*b51+b12*a32*b31*b41+b52*a22-b12*a32*b11*b61+b62*a32+b22*a32*b31*b51+b22*a12*b11*b51-b22*a12*b21*b41-b22*a32*b21*b61+b12*a42+b22*a52+b32*a12*a31*a41-b32*a22*a21*a61+b32*a22*a31*a51+b22*a32*a21*a61-b22*a32*a31*a51-b32*a12*a11*a61+b22*a12*a21*a41-b22*a12*a11*a51; KK[row*30+1*10 + K200 ] = -b22*a32*a11+b12*a22*a31+b22*a12*a31-b12*a32*a21-b32*a12*a21+b32*a22*a11; KK[row*30+1*10 + K110 ] = -2*b12*a12*a31+2*b12*a32*a11-2*b22*a32*a21+2*b22*a22*a31; KK[row*30+1*10 + K020 ] = -b22*a12*a31+b12*a32*a21-b12*a22*a31+b32*a22*a11-b32*a12*a21+b22*a32*a11; KK[row*30+1*10 + K101 ] = -2*b12*a22*a11+2*b12*a12*a21+2*b32*a22*a31-2*b32*a32*a21; KK[row*30+1*10 + K011 ] = 2*b32*a32*a11+2*b22*a12*a21-2*b22*a22*a11-2*b32*a12*a31; KK[row*30+1*10 + K002 ] = -b22*a32*a11+b22*a12*a31+b12*a32*a21+b32*a12*a21-b32*a22*a11-b12*a22*a31; KK[row*30+1*10 + K100 ] = 2*b32*a12*a31-2*b32*a32*a11+2*b22*a12*a21-2*b22*a22*a11; KK[row*30+1*10 + K010 ] = 2*b32*a22*a31-2*b32*a32*a21-2*b12*a12*a21+2*b12*a22*a11; KK[row*30+1*10 + K001 ] = 2*b22*a32*a21-2*b12*a12*a31-2*b22*a22*a31+2*b12*a32*a11; KK[row*30+1*10 + K000 ] = -b32*a22*a11+b22*a32*a11+b32*a12*a21-b22*a12*a31-b12*a32*a21+b12*a22*a31; KK[row*30+2*10 + K200 ] = b22*a12*b31+b22*a32*b11-b12*a32*b21+b12*a22*b31-b32*a22*b11-b32*a12*b21; KK[row*30+2*10 + K110 ] = -2*b12*a12*b31+2*b22*a22*b31-2*b32*a22*b21+2*b32*a12*b11; KK[row*30+2*10 + K020 ] = -b22*a12*b31+b22*a32*b11-b12*a22*b31-b12*a32*b21+b32*a12*b21+b32*a22*b11; KK[row*30+2*10 + K101 ] = 2*b12*a12*b21+2*b22*a32*b31-2*b22*a12*b11-2*b32*a32*b21; KK[row*30+2*10 + K011 ] = -2*b12*a32*b31+2*b12*a22*b21-2*b22*a22*b11+2*b32*a32*b11; KK[row*30+2*10 + K002 ] = -b22*a32*b11-b22*a12*b31+b32*a12*b21+b12*a22*b31+b12*a32*b21-b32*a22*b11; KK[row*30+2*10 + K100 ] = -2*b12*a32*b31-2*b12*a22*b21+2*b22*a22*b11+2*b32*a32*b11; KK[row*30+2*10 + K010 ] = 2*b32*a32*b21-2*b22*a32*b31+2*b12*a12*b21-2*b22*a12*b11; KK[row*30+2*10 + K001 ] = -2*b32*a12*b11+2*b12*a12*b31+2*b22*a22*b31-2*b32*a22*b21; KK[row*30+2*10 + K000 ] = b12*a32*b21+b22*a12*b31-b12*a22*b31-b22*a32*b11+b32*a22*b11-b32*a12*b21; } } #include "solve26AFILE.c" } } #endif void build_last_gb( double *restrow, double *U, int N){ int i,j; double c017, c008, t; /*Här tänkte jag använda V1 och It för att klara en multiplikation! */ for( j=0; j< 84; j++){ restrow[ V3[j] ] += U[32+N*It[j]]; restrow[ V2[j] ] -= U[34+N*It[j]]; } c017 = -U[34+N*56]; c008 = U[32+N*56]; c008 -= c017*U[34+N*56]; for( j=0; j< 84; j++){ restrow[V3[j]] -= c017* U[34 + N*It[j]]; } for( i=0; i< 56; i++){ t = restrow[i]; for( j=0; j< 120; j++) restrow[j] -= t*U[i+N*j]; } for( j=0; j< 120; j++) restrow[j] = restrow[j]/c008; } void build_first_row_of_M( double *restrow, double *U ,double *last_gb, int N){ double v38coeff; double *temp,t; int i,j; v38coeff = - U [33 +N* 56 ]; for( j=0; j< 84; j++){ restrow[ V3[j] ] = -U[33+N*It[j]]; } for( i=0; i< 56; i++){ t = restrow[i]; for( j=0; j< 120; j++) restrow[j] -= t*U[i+N*j]; } for( j=0; j< 120; j++){ restrow[ j ] -= v38coeff*last_gb[j]; } }