subview_cube_meat.hpp 63 KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192939495969798991001011021031041051061071081091101111121131141151161171181191201211221231241251261271281291301311321331341351361371381391401411421431441451461471481491501511521531541551561571581591601611621631641651661671681691701711721731741751761771781791801811821831841851861871881891901911921931941951961971981992002012022032042052062072082092102112122132142152162172182192202212222232242252262272282292302312322332342352362372382392402412422432442452462472482492502512522532542552562572582592602612622632642652662672682692702712722732742752762772782792802812822832842852862872882892902912922932942952962972982993003013023033043053063073083093103113123133143153163173183193203213223233243253263273283293303313323333343353363373383393403413423433443453463473483493503513523533543553563573583593603613623633643653663673683693703713723733743753763773783793803813823833843853863873883893903913923933943953963973983994004014024034044054064074084094104114124134144154164174184194204214224234244254264274284294304314324334344354364374384394404414424434444454464474484494504514524534544554564574584594604614624634644654664674684694704714724734744754764774784794804814824834844854864874884894904914924934944954964974984995005015025035045055065075085095105115125135145155165175185195205215225235245255265275285295305315325335345355365375385395405415425435445455465475485495505515525535545555565575585595605615625635645655665675685695705715725735745755765775785795805815825835845855865875885895905915925935945955965975985996006016026036046056066076086096106116126136146156166176186196206216226236246256266276286296306316326336346356366376386396406416426436446456466476486496506516526536546556566576586596606616626636646656666676686696706716726736746756766776786796806816826836846856866876886896906916926936946956966976986997007017027037047057067077087097107117127137147157167177187197207217227237247257267277287297307317327337347357367377387397407417427437447457467477487497507517527537547557567577587597607617627637647657667677687697707717727737747757767777787797807817827837847857867877887897907917927937947957967977987998008018028038048058068078088098108118128138148158168178188198208218228238248258268278288298308318328338348358368378388398408418428438448458468478488498508518528538548558568578588598608618628638648658668678688698708718728738748758768778788798808818828838848858868878888898908918928938948958968978988999009019029039049059069079089099109119129139149159169179189199209219229239249259269279289299309319329339349359369379389399409419429439449459469479489499509519529539549559569579589599609619629639649659669679689699709719729739749759769779789799809819829839849859869879889899909919929939949959969979989991000100110021003100410051006100710081009101010111012101310141015101610171018101910201021102210231024102510261027102810291030103110321033103410351036103710381039104010411042104310441045104610471048104910501051105210531054105510561057105810591060106110621063106410651066106710681069107010711072107310741075107610771078107910801081108210831084108510861087108810891090109110921093109410951096109710981099110011011102110311041105110611071108110911101111111211131114111511161117111811191120112111221123112411251126112711281129113011311132113311341135113611371138113911401141114211431144114511461147114811491150115111521153115411551156115711581159116011611162116311641165116611671168116911701171117211731174117511761177117811791180118111821183118411851186118711881189119011911192119311941195119611971198119912001201120212031204120512061207120812091210121112121213121412151216121712181219122012211222122312241225122612271228122912301231123212331234123512361237123812391240124112421243124412451246124712481249125012511252125312541255125612571258125912601261126212631264126512661267126812691270127112721273127412751276127712781279128012811282128312841285128612871288128912901291129212931294129512961297129812991300130113021303130413051306130713081309131013111312131313141315131613171318131913201321132213231324132513261327132813291330133113321333133413351336133713381339134013411342134313441345134613471348134913501351135213531354135513561357135813591360136113621363136413651366136713681369137013711372137313741375137613771378137913801381138213831384138513861387138813891390139113921393139413951396139713981399140014011402140314041405140614071408140914101411141214131414141514161417141814191420142114221423142414251426142714281429143014311432143314341435143614371438143914401441144214431444144514461447144814491450145114521453145414551456145714581459146014611462146314641465146614671468146914701471147214731474147514761477147814791480148114821483148414851486148714881489149014911492149314941495149614971498149915001501150215031504150515061507150815091510151115121513151415151516151715181519152015211522152315241525152615271528152915301531153215331534153515361537153815391540154115421543154415451546154715481549155015511552155315541555155615571558155915601561156215631564156515661567156815691570157115721573157415751576157715781579158015811582158315841585158615871588158915901591159215931594159515961597159815991600160116021603160416051606160716081609161016111612161316141615161616171618161916201621162216231624162516261627162816291630163116321633163416351636163716381639164016411642164316441645164616471648164916501651165216531654165516561657165816591660166116621663166416651666166716681669167016711672167316741675167616771678167916801681168216831684168516861687168816891690169116921693169416951696169716981699170017011702170317041705170617071708170917101711171217131714171517161717171817191720172117221723172417251726172717281729173017311732173317341735173617371738173917401741174217431744174517461747174817491750175117521753175417551756175717581759176017611762176317641765176617671768176917701771177217731774177517761777177817791780178117821783178417851786178717881789179017911792179317941795179617971798179918001801180218031804180518061807180818091810181118121813181418151816181718181819182018211822182318241825182618271828182918301831183218331834183518361837183818391840184118421843184418451846184718481849185018511852185318541855185618571858185918601861186218631864186518661867186818691870187118721873187418751876187718781879188018811882188318841885188618871888188918901891189218931894189518961897189818991900190119021903190419051906190719081909191019111912191319141915191619171918191919201921192219231924192519261927192819291930193119321933193419351936193719381939194019411942194319441945194619471948194919501951195219531954195519561957195819591960196119621963196419651966196719681969197019711972197319741975197619771978197919801981198219831984198519861987198819891990199119921993199419951996199719981999200020012002200320042005200620072008200920102011201220132014201520162017201820192020202120222023202420252026202720282029203020312032203320342035203620372038203920402041204220432044204520462047204820492050205120522053205420552056205720582059206020612062206320642065206620672068206920702071207220732074207520762077207820792080208120822083208420852086208720882089209020912092209320942095209620972098209921002101210221032104210521062107210821092110211121122113211421152116211721182119212021212122212321242125212621272128212921302131213221332134213521362137213821392140214121422143214421452146214721482149215021512152215321542155215621572158215921602161216221632164216521662167216821692170217121722173217421752176217721782179218021812182218321842185218621872188218921902191219221932194219521962197219821992200220122022203220422052206220722082209221022112212221322142215221622172218221922202221222222232224222522262227222822292230223122322233223422352236223722382239224022412242224322442245224622472248224922502251225222532254225522562257225822592260226122622263226422652266226722682269227022712272227322742275227622772278227922802281228222832284228522862287228822892290229122922293229422952296229722982299230023012302230323042305230623072308230923102311231223132314231523162317231823192320232123222323232423252326232723282329233023312332233323342335233623372338233923402341234223432344234523462347234823492350235123522353235423552356235723582359236023612362236323642365236623672368236923702371237223732374237523762377237823792380238123822383238423852386238723882389239023912392239323942395239623972398239924002401240224032404240524062407240824092410241124122413241424152416241724182419242024212422242324242425242624272428242924302431243224332434243524362437243824392440244124422443244424452446244724482449245024512452245324542455245624572458245924602461246224632464246524662467246824692470247124722473247424752476247724782479248024812482248324842485248624872488248924902491249224932494249524962497249824992500250125022503250425052506250725082509251025112512251325142515251625172518251925202521252225232524252525262527252825292530253125322533253425352536253725382539254025412542254325442545254625472548254925502551255225532554255525562557255825592560256125622563256425652566256725682569257025712572257325742575257625772578257925802581258225832584258525862587258825892590259125922593259425952596259725982599260026012602260326042605260626072608260926102611261226132614261526162617261826192620262126222623262426252626262726282629263026312632263326342635263626372638263926402641264226432644264526462647264826492650265126522653265426552656265726582659266026612662266326642665266626672668266926702671267226732674267526762677267826792680268126822683268426852686268726882689269026912692269326942695269626972698269927002701
  1. // Copyright 2008-2016 Conrad Sanderson (http://conradsanderson.id.au)
  2. // Copyright 2008-2016 National ICT Australia (NICTA)
  3. //
  4. // Licensed under the Apache License, Version 2.0 (the "License");
  5. // you may not use this file except in compliance with the License.
  6. // You may obtain a copy of the License at
  7. // http://www.apache.org/licenses/LICENSE-2.0
  8. //
  9. // Unless required by applicable law or agreed to in writing, software
  10. // distributed under the License is distributed on an "AS IS" BASIS,
  11. // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  12. // See the License for the specific language governing permissions and
  13. // limitations under the License.
  14. // ------------------------------------------------------------------------
  15. //! \addtogroup subview_cube
  16. //! @{
  17. template<typename eT>
  18. inline
  19. subview_cube<eT>::~subview_cube()
  20. {
  21. arma_extra_debug_sigprint();
  22. }
  23. template<typename eT>
  24. arma_inline
  25. subview_cube<eT>::subview_cube
  26. (
  27. const Cube<eT>& in_m,
  28. const uword in_row1,
  29. const uword in_col1,
  30. const uword in_slice1,
  31. const uword in_n_rows,
  32. const uword in_n_cols,
  33. const uword in_n_slices
  34. )
  35. : m (in_m)
  36. , aux_row1 (in_row1)
  37. , aux_col1 (in_col1)
  38. , aux_slice1 (in_slice1)
  39. , n_rows (in_n_rows)
  40. , n_cols (in_n_cols)
  41. , n_elem_slice(in_n_rows * in_n_cols)
  42. , n_slices (in_n_slices)
  43. , n_elem (n_elem_slice * in_n_slices)
  44. {
  45. arma_extra_debug_sigprint();
  46. }
  47. template<typename eT>
  48. inline
  49. void
  50. subview_cube<eT>::operator= (const eT val)
  51. {
  52. arma_extra_debug_sigprint();
  53. if(n_elem != 1)
  54. {
  55. arma_debug_assert_same_size(n_rows, n_cols, n_slices, 1, 1, 1, "copy into subcube");
  56. }
  57. Cube<eT>& Q = const_cast< Cube<eT>& >(m);
  58. Q.at(aux_row1, aux_col1, aux_slice1) = val;
  59. }
  60. template<typename eT>
  61. inline
  62. void
  63. subview_cube<eT>::operator+= (const eT val)
  64. {
  65. arma_extra_debug_sigprint();
  66. const uword local_n_rows = n_rows;
  67. const uword local_n_cols = n_cols;
  68. const uword local_n_slices = n_slices;
  69. for(uword slice = 0; slice < local_n_slices; ++slice)
  70. {
  71. for(uword col = 0; col < local_n_cols; ++col)
  72. {
  73. arrayops::inplace_plus( slice_colptr(slice,col), val, local_n_rows );
  74. }
  75. }
  76. }
  77. template<typename eT>
  78. inline
  79. void
  80. subview_cube<eT>::operator-= (const eT val)
  81. {
  82. arma_extra_debug_sigprint();
  83. const uword local_n_rows = n_rows;
  84. const uword local_n_cols = n_cols;
  85. const uword local_n_slices = n_slices;
  86. for(uword slice = 0; slice < local_n_slices; ++slice)
  87. {
  88. for(uword col = 0; col < local_n_cols; ++col)
  89. {
  90. arrayops::inplace_minus( slice_colptr(slice,col), val, local_n_rows );
  91. }
  92. }
  93. }
  94. template<typename eT>
  95. inline
  96. void
  97. subview_cube<eT>::operator*= (const eT val)
  98. {
  99. arma_extra_debug_sigprint();
  100. const uword local_n_rows = n_rows;
  101. const uword local_n_cols = n_cols;
  102. const uword local_n_slices = n_slices;
  103. for(uword slice = 0; slice < local_n_slices; ++slice)
  104. {
  105. for(uword col = 0; col < local_n_cols; ++col)
  106. {
  107. arrayops::inplace_mul( slice_colptr(slice,col), val, local_n_rows );
  108. }
  109. }
  110. }
  111. template<typename eT>
  112. inline
  113. void
  114. subview_cube<eT>::operator/= (const eT val)
  115. {
  116. arma_extra_debug_sigprint();
  117. const uword local_n_rows = n_rows;
  118. const uword local_n_cols = n_cols;
  119. const uword local_n_slices = n_slices;
  120. for(uword slice = 0; slice < local_n_slices; ++slice)
  121. {
  122. for(uword col = 0; col < local_n_cols; ++col)
  123. {
  124. arrayops::inplace_div( slice_colptr(slice,col), val, local_n_rows );
  125. }
  126. }
  127. }
  128. template<typename eT>
  129. template<typename T1>
  130. inline
  131. void
  132. subview_cube<eT>::operator= (const BaseCube<eT,T1>& in)
  133. {
  134. arma_extra_debug_sigprint();
  135. const unwrap_cube<T1> tmp(in.get_ref());
  136. const Cube<eT>& x = tmp.M;
  137. subview_cube<eT>& t = *this;
  138. arma_debug_assert_same_size(t, x, "copy into subcube");
  139. const uword t_n_rows = t.n_rows;
  140. const uword t_n_cols = t.n_cols;
  141. const uword t_n_slices = t.n_slices;
  142. for(uword slice = 0; slice < t_n_slices; ++slice)
  143. {
  144. for(uword col = 0; col < t_n_cols; ++col)
  145. {
  146. arrayops::copy( t.slice_colptr(slice,col), x.slice_colptr(slice,col), t_n_rows );
  147. }
  148. }
  149. }
  150. template<typename eT>
  151. template<typename T1>
  152. inline
  153. void
  154. subview_cube<eT>::operator+= (const BaseCube<eT,T1>& in)
  155. {
  156. arma_extra_debug_sigprint();
  157. const unwrap_cube<T1> tmp(in.get_ref());
  158. const Cube<eT>& x = tmp.M;
  159. subview_cube<eT>& t = *this;
  160. arma_debug_assert_same_size(t, x, "addition");
  161. const uword t_n_rows = t.n_rows;
  162. const uword t_n_cols = t.n_cols;
  163. const uword t_n_slices = t.n_slices;
  164. for(uword slice = 0; slice < t_n_slices; ++slice)
  165. {
  166. for(uword col = 0; col < t_n_cols; ++col)
  167. {
  168. arrayops::inplace_plus( t.slice_colptr(slice,col), x.slice_colptr(slice,col), t_n_rows );
  169. }
  170. }
  171. }
  172. template<typename eT>
  173. template<typename T1>
  174. inline
  175. void
  176. subview_cube<eT>::operator-= (const BaseCube<eT,T1>& in)
  177. {
  178. arma_extra_debug_sigprint();
  179. const unwrap_cube<T1> tmp(in.get_ref());
  180. const Cube<eT>& x = tmp.M;
  181. subview_cube<eT>& t = *this;
  182. arma_debug_assert_same_size(t, x, "subtraction");
  183. const uword t_n_rows = t.n_rows;
  184. const uword t_n_cols = t.n_cols;
  185. const uword t_n_slices = t.n_slices;
  186. for(uword slice = 0; slice < t_n_slices; ++slice)
  187. {
  188. for(uword col = 0; col < t_n_cols; ++col)
  189. {
  190. arrayops::inplace_minus( t.slice_colptr(slice,col), x.slice_colptr(slice,col), t_n_rows );
  191. }
  192. }
  193. }
  194. template<typename eT>
  195. template<typename T1>
  196. inline
  197. void
  198. subview_cube<eT>::operator%= (const BaseCube<eT,T1>& in)
  199. {
  200. arma_extra_debug_sigprint();
  201. const unwrap_cube<T1> tmp(in.get_ref());
  202. const Cube<eT>& x = tmp.M;
  203. subview_cube<eT>& t = *this;
  204. arma_debug_assert_same_size(t, x, "element-wise multiplication");
  205. const uword t_n_rows = t.n_rows;
  206. const uword t_n_cols = t.n_cols;
  207. const uword t_n_slices = t.n_slices;
  208. for(uword slice = 0; slice < t_n_slices; ++slice)
  209. {
  210. for(uword col = 0; col < t_n_cols; ++col)
  211. {
  212. arrayops::inplace_mul( t.slice_colptr(slice,col), x.slice_colptr(slice,col), t_n_rows );
  213. }
  214. }
  215. }
  216. template<typename eT>
  217. template<typename T1>
  218. inline
  219. void
  220. subview_cube<eT>::operator/= (const BaseCube<eT,T1>& in)
  221. {
  222. arma_extra_debug_sigprint();
  223. const unwrap_cube<T1> tmp(in.get_ref());
  224. const Cube<eT>& x = tmp.M;
  225. subview_cube<eT>& t = *this;
  226. arma_debug_assert_same_size(t, x, "element-wise division");
  227. const uword t_n_rows = t.n_rows;
  228. const uword t_n_cols = t.n_cols;
  229. const uword t_n_slices = t.n_slices;
  230. for(uword slice = 0; slice < t_n_slices; ++slice)
  231. {
  232. for(uword col = 0; col < t_n_cols; ++col)
  233. {
  234. arrayops::inplace_div( t.slice_colptr(slice,col), x.slice_colptr(slice,col), t_n_rows );
  235. }
  236. }
  237. }
  238. //! x.subcube(...) = y.subcube(...)
  239. template<typename eT>
  240. inline
  241. void
  242. subview_cube<eT>::operator= (const subview_cube<eT>& x)
  243. {
  244. arma_extra_debug_sigprint();
  245. if(check_overlap(x))
  246. {
  247. const Cube<eT> tmp(x);
  248. (*this).operator=(tmp);
  249. return;
  250. }
  251. subview_cube<eT>& t = *this;
  252. arma_debug_assert_same_size(t, x, "copy into subcube");
  253. const uword t_n_rows = t.n_rows;
  254. const uword t_n_cols = t.n_cols;
  255. const uword t_n_slices = t.n_slices;
  256. for(uword slice = 0; slice < t_n_slices; ++slice)
  257. {
  258. for(uword col = 0; col < t_n_cols; ++col)
  259. {
  260. arrayops::copy( t.slice_colptr(slice,col), x.slice_colptr(slice,col), t_n_rows );
  261. }
  262. }
  263. }
  264. template<typename eT>
  265. inline
  266. void
  267. subview_cube<eT>::operator+= (const subview_cube<eT>& x)
  268. {
  269. arma_extra_debug_sigprint();
  270. if(check_overlap(x))
  271. {
  272. const Cube<eT> tmp(x);
  273. (*this).operator+=(tmp);
  274. return;
  275. }
  276. subview_cube<eT>& t = *this;
  277. arma_debug_assert_same_size(t, x, "addition");
  278. const uword t_n_rows = t.n_rows;
  279. const uword t_n_cols = t.n_cols;
  280. const uword t_n_slices = t.n_slices;
  281. for(uword slice = 0; slice < t_n_slices; ++slice)
  282. {
  283. for(uword col = 0; col < t_n_cols; ++col)
  284. {
  285. arrayops::inplace_plus( t.slice_colptr(slice,col), x.slice_colptr(slice,col), t_n_rows );
  286. }
  287. }
  288. }
  289. template<typename eT>
  290. inline
  291. void
  292. subview_cube<eT>::operator-= (const subview_cube<eT>& x)
  293. {
  294. arma_extra_debug_sigprint();
  295. if(check_overlap(x))
  296. {
  297. const Cube<eT> tmp(x);
  298. (*this).operator-=(tmp);
  299. return;
  300. }
  301. subview_cube<eT>& t = *this;
  302. arma_debug_assert_same_size(t, x, "subtraction");
  303. const uword t_n_rows = t.n_rows;
  304. const uword t_n_cols = t.n_cols;
  305. const uword t_n_slices = t.n_slices;
  306. for(uword slice = 0; slice < t_n_slices; ++slice)
  307. {
  308. for(uword col = 0; col < t_n_cols; ++col)
  309. {
  310. arrayops::inplace_minus( t.slice_colptr(slice,col), x.slice_colptr(slice,col), t_n_rows );
  311. }
  312. }
  313. }
  314. template<typename eT>
  315. inline
  316. void
  317. subview_cube<eT>::operator%= (const subview_cube<eT>& x)
  318. {
  319. arma_extra_debug_sigprint();
  320. if(check_overlap(x))
  321. {
  322. const Cube<eT> tmp(x);
  323. (*this).operator%=(tmp);
  324. return;
  325. }
  326. subview_cube<eT>& t = *this;
  327. arma_debug_assert_same_size(t, x, "element-wise multiplication");
  328. const uword t_n_rows = t.n_rows;
  329. const uword t_n_cols = t.n_cols;
  330. const uword t_n_slices = t.n_slices;
  331. for(uword slice = 0; slice < t_n_slices; ++slice)
  332. {
  333. for(uword col = 0; col < t_n_cols; ++col)
  334. {
  335. arrayops::inplace_mul( t.slice_colptr(slice,col), x.slice_colptr(slice,col), t_n_rows );
  336. }
  337. }
  338. }
  339. template<typename eT>
  340. inline
  341. void
  342. subview_cube<eT>::operator/= (const subview_cube<eT>& x)
  343. {
  344. arma_extra_debug_sigprint();
  345. if(check_overlap(x))
  346. {
  347. const Cube<eT> tmp(x);
  348. (*this).operator/=(tmp);
  349. return;
  350. }
  351. subview_cube<eT>& t = *this;
  352. arma_debug_assert_same_size(t, x, "element-wise division");
  353. const uword t_n_rows = t.n_rows;
  354. const uword t_n_cols = t.n_cols;
  355. const uword t_n_slices = t.n_slices;
  356. for(uword slice = 0; slice < t_n_slices; ++slice)
  357. {
  358. for(uword col = 0; col < t_n_cols; ++col)
  359. {
  360. arrayops::inplace_div( t.slice_colptr(slice,col), x.slice_colptr(slice,col), t_n_rows );
  361. }
  362. }
  363. }
  364. template<typename eT>
  365. template<typename T1>
  366. inline
  367. void
  368. subview_cube<eT>::operator= (const Base<eT,T1>& in)
  369. {
  370. arma_extra_debug_sigprint();
  371. const quasi_unwrap<T1> tmp(in.get_ref());
  372. const Mat<eT>& x = tmp.M;
  373. subview_cube<eT>& t = *this;
  374. const uword t_n_rows = t.n_rows;
  375. const uword t_n_cols = t.n_cols;
  376. const uword t_n_slices = t.n_slices;
  377. const uword x_n_rows = x.n_rows;
  378. const uword x_n_cols = x.n_cols;
  379. if( ((x_n_rows == 1) || (x_n_cols == 1)) && (t_n_rows == 1) && (t_n_cols == 1) && (x.n_elem == t_n_slices) )
  380. {
  381. Cube<eT>& Q = const_cast< Cube<eT>& >(t.m);
  382. const uword t_aux_row1 = t.aux_row1;
  383. const uword t_aux_col1 = t.aux_col1;
  384. const uword t_aux_slice1 = t.aux_slice1;
  385. const eT* x_mem = x.memptr();
  386. uword i,j;
  387. for(i=0, j=1; j < t_n_slices; i+=2, j+=2)
  388. {
  389. const eT tmp_i = x_mem[i];
  390. const eT tmp_j = x_mem[j];
  391. Q.at(t_aux_row1, t_aux_col1, t_aux_slice1 + i) = tmp_i;
  392. Q.at(t_aux_row1, t_aux_col1, t_aux_slice1 + j) = tmp_j;
  393. }
  394. if(i < t_n_slices)
  395. {
  396. Q.at(t_aux_row1, t_aux_col1, t_aux_slice1 + i) = x_mem[i];
  397. }
  398. }
  399. else
  400. if( (t_n_rows == x_n_rows) && (t_n_cols == x_n_cols) && (t_n_slices == 1) )
  401. {
  402. // interpret the matrix as a cube with one slice
  403. for(uword col = 0; col < t_n_cols; ++col)
  404. {
  405. arrayops::copy( t.slice_colptr(0, col), x.colptr(col), t_n_rows );
  406. }
  407. }
  408. else
  409. if( (t_n_rows == x_n_rows) && (t_n_cols == 1) && (t_n_slices == x_n_cols) )
  410. {
  411. for(uword i=0; i < t_n_slices; ++i)
  412. {
  413. arrayops::copy( t.slice_colptr(i, 0), x.colptr(i), t_n_rows );
  414. }
  415. }
  416. else
  417. if( (t_n_rows == 1) && (t_n_cols == x_n_rows) && (t_n_slices == x_n_cols) )
  418. {
  419. Cube<eT>& Q = const_cast< Cube<eT>& >(t.m);
  420. const uword t_aux_row1 = t.aux_row1;
  421. const uword t_aux_col1 = t.aux_col1;
  422. const uword t_aux_slice1 = t.aux_slice1;
  423. for(uword slice=0; slice < t_n_slices; ++slice)
  424. {
  425. const uword mod_slice = t_aux_slice1 + slice;
  426. const eT* x_colptr = x.colptr(slice);
  427. uword i,j;
  428. for(i=0, j=1; j < t_n_cols; i+=2, j+=2)
  429. {
  430. const eT tmp_i = x_colptr[i];
  431. const eT tmp_j = x_colptr[j];
  432. Q.at(t_aux_row1, t_aux_col1 + i, mod_slice) = tmp_i;
  433. Q.at(t_aux_row1, t_aux_col1 + j, mod_slice) = tmp_j;
  434. }
  435. if(i < t_n_cols)
  436. {
  437. Q.at(t_aux_row1, t_aux_col1 + i, mod_slice) = x_colptr[i];
  438. }
  439. }
  440. }
  441. else
  442. {
  443. if(arma_config::debug)
  444. {
  445. arma_stop_logic_error( arma_incompat_size_string(t, x, "copy into subcube") );
  446. }
  447. }
  448. }
  449. template<typename eT>
  450. template<typename T1>
  451. inline
  452. void
  453. subview_cube<eT>::operator+= (const Base<eT,T1>& in)
  454. {
  455. arma_extra_debug_sigprint();
  456. const quasi_unwrap<T1> tmp(in.get_ref());
  457. const Mat<eT>& x = tmp.M;
  458. subview_cube<eT>& t = *this;
  459. const uword t_n_rows = t.n_rows;
  460. const uword t_n_cols = t.n_cols;
  461. const uword t_n_slices = t.n_slices;
  462. const uword x_n_rows = x.n_rows;
  463. const uword x_n_cols = x.n_cols;
  464. if( ((x_n_rows == 1) || (x_n_cols == 1)) && (t_n_rows == 1) && (t_n_cols == 1) && (x.n_elem == t_n_slices) )
  465. {
  466. Cube<eT>& Q = const_cast< Cube<eT>& >(t.m);
  467. const uword t_aux_row1 = t.aux_row1;
  468. const uword t_aux_col1 = t.aux_col1;
  469. const uword t_aux_slice1 = t.aux_slice1;
  470. const eT* x_mem = x.memptr();
  471. uword i,j;
  472. for(i=0, j=1; j < t_n_slices; i+=2, j+=2)
  473. {
  474. const eT tmp_i = x_mem[i];
  475. const eT tmp_j = x_mem[j];
  476. Q.at(t_aux_row1, t_aux_col1, t_aux_slice1 + i) += tmp_i;
  477. Q.at(t_aux_row1, t_aux_col1, t_aux_slice1 + j) += tmp_j;
  478. }
  479. if(i < t_n_slices)
  480. {
  481. Q.at(t_aux_row1, t_aux_col1, t_aux_slice1 + i) += x_mem[i];
  482. }
  483. }
  484. else
  485. if( (t_n_rows == x_n_rows) && (t_n_cols == x_n_cols) && (t_n_slices == 1) )
  486. {
  487. for(uword col = 0; col < t_n_cols; ++col)
  488. {
  489. arrayops::inplace_plus( t.slice_colptr(0, col), x.colptr(col), t_n_rows );
  490. }
  491. }
  492. else
  493. if( (t_n_rows == x_n_rows) && (t_n_cols == 1) && (t_n_slices == x_n_cols) )
  494. {
  495. for(uword i=0; i < t_n_slices; ++i)
  496. {
  497. arrayops::inplace_plus( t.slice_colptr(i, 0), x.colptr(i), t_n_rows );
  498. }
  499. }
  500. else
  501. if( (t_n_rows == 1) && (t_n_cols == x_n_rows) && (t_n_slices == x_n_cols) )
  502. {
  503. Cube<eT>& Q = const_cast< Cube<eT>& >(t.m);
  504. const uword t_aux_row1 = t.aux_row1;
  505. const uword t_aux_col1 = t.aux_col1;
  506. const uword t_aux_slice1 = t.aux_slice1;
  507. for(uword slice=0; slice < t_n_slices; ++slice)
  508. {
  509. const uword mod_slice = t_aux_slice1 + slice;
  510. const eT* x_colptr = x.colptr(slice);
  511. uword i,j;
  512. for(i=0, j=1; j < t_n_cols; i+=2, j+=2)
  513. {
  514. const eT tmp_i = x_colptr[i];
  515. const eT tmp_j = x_colptr[j];
  516. Q.at(t_aux_row1, t_aux_col1 + i, mod_slice) += tmp_i;
  517. Q.at(t_aux_row1, t_aux_col1 + j, mod_slice) += tmp_j;
  518. }
  519. if(i < t_n_cols)
  520. {
  521. Q.at(t_aux_row1, t_aux_col1 + i, mod_slice) += x_colptr[i];
  522. }
  523. }
  524. }
  525. else
  526. {
  527. if(arma_config::debug)
  528. {
  529. arma_stop_logic_error( arma_incompat_size_string(t, x, "addition") );
  530. }
  531. }
  532. }
  533. template<typename eT>
  534. template<typename T1>
  535. inline
  536. void
  537. subview_cube<eT>::operator-= (const Base<eT,T1>& in)
  538. {
  539. arma_extra_debug_sigprint();
  540. const quasi_unwrap<T1> tmp(in.get_ref());
  541. const Mat<eT>& x = tmp.M;
  542. subview_cube<eT>& t = *this;
  543. const uword t_n_rows = t.n_rows;
  544. const uword t_n_cols = t.n_cols;
  545. const uword t_n_slices = t.n_slices;
  546. const uword x_n_rows = x.n_rows;
  547. const uword x_n_cols = x.n_cols;
  548. if( ((x_n_rows == 1) || (x_n_cols == 1)) && (t_n_rows == 1) && (t_n_cols == 1) && (x.n_elem == t_n_slices) )
  549. {
  550. Cube<eT>& Q = const_cast< Cube<eT>& >(t.m);
  551. const uword t_aux_row1 = t.aux_row1;
  552. const uword t_aux_col1 = t.aux_col1;
  553. const uword t_aux_slice1 = t.aux_slice1;
  554. const eT* x_mem = x.memptr();
  555. uword i,j;
  556. for(i=0, j=1; j < t_n_slices; i+=2, j+=2)
  557. {
  558. const eT tmp_i = x_mem[i];
  559. const eT tmp_j = x_mem[j];
  560. Q.at(t_aux_row1, t_aux_col1, t_aux_slice1 + i) -= tmp_i;
  561. Q.at(t_aux_row1, t_aux_col1, t_aux_slice1 + j) -= tmp_j;
  562. }
  563. if(i < t_n_slices)
  564. {
  565. Q.at(t_aux_row1, t_aux_col1, t_aux_slice1 + i) -= x_mem[i];
  566. }
  567. }
  568. else
  569. if( (t_n_rows == x_n_rows) && (t_n_cols == x_n_cols) && (t_n_slices == 1) )
  570. {
  571. for(uword col = 0; col < t_n_cols; ++col)
  572. {
  573. arrayops::inplace_minus( t.slice_colptr(0, col), x.colptr(col), t_n_rows );
  574. }
  575. }
  576. else
  577. if( (t_n_rows == x_n_rows) && (t_n_cols == 1) && (t_n_slices == x_n_cols) )
  578. {
  579. for(uword i=0; i < t_n_slices; ++i)
  580. {
  581. arrayops::inplace_minus( t.slice_colptr(i, 0), x.colptr(i), t_n_rows );
  582. }
  583. }
  584. else
  585. if( (t_n_rows == 1) && (t_n_cols == x_n_rows) && (t_n_slices == x_n_cols) )
  586. {
  587. Cube<eT>& Q = const_cast< Cube<eT>& >(t.m);
  588. const uword t_aux_row1 = t.aux_row1;
  589. const uword t_aux_col1 = t.aux_col1;
  590. const uword t_aux_slice1 = t.aux_slice1;
  591. for(uword slice=0; slice < t_n_slices; ++slice)
  592. {
  593. const uword mod_slice = t_aux_slice1 + slice;
  594. const eT* x_colptr = x.colptr(slice);
  595. uword i,j;
  596. for(i=0, j=1; j < t_n_cols; i+=2, j+=2)
  597. {
  598. const eT tmp_i = x_colptr[i];
  599. const eT tmp_j = x_colptr[j];
  600. Q.at(t_aux_row1, t_aux_col1 + i, mod_slice) -= tmp_i;
  601. Q.at(t_aux_row1, t_aux_col1 + j, mod_slice) -= tmp_j;
  602. }
  603. if(i < t_n_cols)
  604. {
  605. Q.at(t_aux_row1, t_aux_col1 + i, mod_slice) -= x_colptr[i];
  606. }
  607. }
  608. }
  609. else
  610. {
  611. if(arma_config::debug)
  612. {
  613. arma_stop_logic_error( arma_incompat_size_string(t, x, "subtraction") );
  614. }
  615. }
  616. }
  617. template<typename eT>
  618. template<typename T1>
  619. inline
  620. void
  621. subview_cube<eT>::operator%= (const Base<eT,T1>& in)
  622. {
  623. arma_extra_debug_sigprint();
  624. const quasi_unwrap<T1> tmp(in.get_ref());
  625. const Mat<eT>& x = tmp.M;
  626. subview_cube<eT>& t = *this;
  627. const uword t_n_rows = t.n_rows;
  628. const uword t_n_cols = t.n_cols;
  629. const uword t_n_slices = t.n_slices;
  630. const uword x_n_rows = x.n_rows;
  631. const uword x_n_cols = x.n_cols;
  632. if( ((x_n_rows == 1) || (x_n_cols == 1)) && (t_n_rows == 1) && (t_n_cols == 1) && (x.n_elem == t_n_slices) )
  633. {
  634. Cube<eT>& Q = const_cast< Cube<eT>& >(t.m);
  635. const uword t_aux_row1 = t.aux_row1;
  636. const uword t_aux_col1 = t.aux_col1;
  637. const uword t_aux_slice1 = t.aux_slice1;
  638. const eT* x_mem = x.memptr();
  639. uword i,j;
  640. for(i=0, j=1; j < t_n_slices; i+=2, j+=2)
  641. {
  642. const eT tmp_i = x_mem[i];
  643. const eT tmp_j = x_mem[j];
  644. Q.at(t_aux_row1, t_aux_col1, t_aux_slice1 + i) *= tmp_i;
  645. Q.at(t_aux_row1, t_aux_col1, t_aux_slice1 + j) *= tmp_j;
  646. }
  647. if(i < t_n_slices)
  648. {
  649. Q.at(t_aux_row1, t_aux_col1, t_aux_slice1 + i) *= x_mem[i];
  650. }
  651. }
  652. else
  653. if( (t_n_rows == x_n_rows) && (t_n_cols == x_n_cols) && (t_n_slices == 1) )
  654. {
  655. for(uword col = 0; col < t_n_cols; ++col)
  656. {
  657. arrayops::inplace_mul( t.slice_colptr(0, col), x.colptr(col), t_n_rows );
  658. }
  659. }
  660. else
  661. if( (t_n_rows == x_n_rows) && (t_n_cols == 1) && (t_n_slices == x_n_cols) )
  662. {
  663. for(uword i=0; i < t_n_slices; ++i)
  664. {
  665. arrayops::inplace_mul( t.slice_colptr(i, 0), x.colptr(i), t_n_rows );
  666. }
  667. }
  668. else
  669. if( (t_n_rows == 1) && (t_n_cols == x_n_rows) && (t_n_slices == x_n_cols) )
  670. {
  671. Cube<eT>& Q = const_cast< Cube<eT>& >(t.m);
  672. const uword t_aux_row1 = t.aux_row1;
  673. const uword t_aux_col1 = t.aux_col1;
  674. const uword t_aux_slice1 = t.aux_slice1;
  675. for(uword slice=0; slice < t_n_slices; ++slice)
  676. {
  677. const uword mod_slice = t_aux_slice1 + slice;
  678. const eT* x_colptr = x.colptr(slice);
  679. uword i,j;
  680. for(i=0, j=1; j < t_n_cols; i+=2, j+=2)
  681. {
  682. const eT tmp_i = x_colptr[i];
  683. const eT tmp_j = x_colptr[j];
  684. Q.at(t_aux_row1, t_aux_col1 + i, mod_slice) *= tmp_i;
  685. Q.at(t_aux_row1, t_aux_col1 + j, mod_slice) *= tmp_j;
  686. }
  687. if(i < t_n_cols)
  688. {
  689. Q.at(t_aux_row1, t_aux_col1 + i, mod_slice) *= x_colptr[i];
  690. }
  691. }
  692. }
  693. else
  694. {
  695. if(arma_config::debug)
  696. {
  697. arma_stop_logic_error( arma_incompat_size_string(t, x, "element-wise multiplication") );
  698. }
  699. }
  700. }
  701. template<typename eT>
  702. template<typename T1>
  703. inline
  704. void
  705. subview_cube<eT>::operator/= (const Base<eT,T1>& in)
  706. {
  707. arma_extra_debug_sigprint();
  708. const quasi_unwrap<T1> tmp(in.get_ref());
  709. const Mat<eT>& x = tmp.M;
  710. subview_cube<eT>& t = *this;
  711. const uword t_n_rows = t.n_rows;
  712. const uword t_n_cols = t.n_cols;
  713. const uword t_n_slices = t.n_slices;
  714. const uword x_n_rows = x.n_rows;
  715. const uword x_n_cols = x.n_cols;
  716. if( ((x_n_rows == 1) || (x_n_cols == 1)) && (t_n_rows == 1) && (t_n_cols == 1) && (x.n_elem == t_n_slices) )
  717. {
  718. Cube<eT>& Q = const_cast< Cube<eT>& >(t.m);
  719. const uword t_aux_row1 = t.aux_row1;
  720. const uword t_aux_col1 = t.aux_col1;
  721. const uword t_aux_slice1 = t.aux_slice1;
  722. const eT* x_mem = x.memptr();
  723. uword i,j;
  724. for(i=0, j=1; j < t_n_slices; i+=2, j+=2)
  725. {
  726. const eT tmp_i = x_mem[i];
  727. const eT tmp_j = x_mem[j];
  728. Q.at(t_aux_row1, t_aux_col1, t_aux_slice1 + i) /= tmp_i;
  729. Q.at(t_aux_row1, t_aux_col1, t_aux_slice1 + j) /= tmp_j;
  730. }
  731. if(i < t_n_slices)
  732. {
  733. Q.at(t_aux_row1, t_aux_col1, t_aux_slice1 + i) /= x_mem[i];
  734. }
  735. }
  736. else
  737. if( (t_n_rows == x_n_rows) && (t_n_cols == x_n_cols) && (t_n_slices == 1) )
  738. {
  739. for(uword col = 0; col < t_n_cols; ++col)
  740. {
  741. arrayops::inplace_div( t.slice_colptr(0, col), x.colptr(col), t_n_rows );
  742. }
  743. }
  744. else
  745. if( (t_n_rows == x_n_rows) && (t_n_cols == 1) && (t_n_slices == x_n_cols) )
  746. {
  747. for(uword i=0; i < t_n_slices; ++i)
  748. {
  749. arrayops::inplace_div( t.slice_colptr(i, 0), x.colptr(i), t_n_rows );
  750. }
  751. }
  752. else
  753. if( (t_n_rows == 1) && (t_n_cols == x_n_rows) && (t_n_slices == x_n_cols) )
  754. {
  755. Cube<eT>& Q = const_cast< Cube<eT>& >(t.m);
  756. const uword t_aux_row1 = t.aux_row1;
  757. const uword t_aux_col1 = t.aux_col1;
  758. const uword t_aux_slice1 = t.aux_slice1;
  759. for(uword slice=0; slice < t_n_slices; ++slice)
  760. {
  761. const uword mod_slice = t_aux_slice1 + slice;
  762. const eT* x_colptr = x.colptr(slice);
  763. uword i,j;
  764. for(i=0, j=1; j < t_n_cols; i+=2, j+=2)
  765. {
  766. const eT tmp_i = x_colptr[i];
  767. const eT tmp_j = x_colptr[j];
  768. Q.at(t_aux_row1, t_aux_col1 + i, mod_slice) /= tmp_i;
  769. Q.at(t_aux_row1, t_aux_col1 + j, mod_slice) /= tmp_j;
  770. }
  771. if(i < t_n_cols)
  772. {
  773. Q.at(t_aux_row1, t_aux_col1 + i, mod_slice) /= x_colptr[i];
  774. }
  775. }
  776. }
  777. else
  778. {
  779. if(arma_config::debug)
  780. {
  781. arma_stop_logic_error( arma_incompat_size_string(t, x, "element-wise division") );
  782. }
  783. }
  784. }
  785. template<typename eT>
  786. template<typename gen_type>
  787. inline
  788. void
  789. subview_cube<eT>::operator= (const GenCube<eT,gen_type>& in)
  790. {
  791. arma_extra_debug_sigprint();
  792. arma_debug_assert_same_size(n_rows, n_cols, n_slices, in.n_rows, in.n_cols, in.n_slices, "copy into subcube");
  793. in.apply(*this);
  794. }
  795. //! apply a functor to each element
  796. template<typename eT>
  797. template<typename functor>
  798. inline
  799. void
  800. subview_cube<eT>::for_each(functor F)
  801. {
  802. arma_extra_debug_sigprint();
  803. Cube<eT>& Q = const_cast< Cube<eT>& >(m);
  804. const uword start_col = aux_col1;
  805. const uword start_row = aux_row1;
  806. const uword start_slice = aux_slice1;
  807. const uword end_col_plus1 = start_col + n_cols;
  808. const uword end_row_plus1 = start_row + n_rows;
  809. const uword end_slice_plus1 = start_slice + n_slices;
  810. for(uword uslice = start_slice; uslice < end_slice_plus1; ++uslice)
  811. for(uword ucol = start_col; ucol < end_col_plus1; ++ucol )
  812. for(uword urow = start_row; urow < end_row_plus1; ++urow )
  813. {
  814. F( Q.at(urow, ucol, uslice) );
  815. }
  816. }
  817. template<typename eT>
  818. template<typename functor>
  819. inline
  820. void
  821. subview_cube<eT>::for_each(functor F) const
  822. {
  823. arma_extra_debug_sigprint();
  824. const Cube<eT>& Q = m;
  825. const uword start_col = aux_col1;
  826. const uword start_row = aux_row1;
  827. const uword start_slice = aux_slice1;
  828. const uword end_col_plus1 = start_col + n_cols;
  829. const uword end_row_plus1 = start_row + n_rows;
  830. const uword end_slice_plus1 = start_slice + n_slices;
  831. for(uword uslice = start_slice; uslice < end_slice_plus1; ++uslice)
  832. for(uword ucol = start_col; ucol < end_col_plus1; ++ucol )
  833. for(uword urow = start_row; urow < end_row_plus1; ++urow )
  834. {
  835. F( Q.at(urow, ucol, uslice) );
  836. }
  837. }
  838. //! transform each element in the subview using a functor
  839. template<typename eT>
  840. template<typename functor>
  841. inline
  842. void
  843. subview_cube<eT>::transform(functor F)
  844. {
  845. arma_extra_debug_sigprint();
  846. Cube<eT>& Q = const_cast< Cube<eT>& >(m);
  847. const uword start_col = aux_col1;
  848. const uword start_row = aux_row1;
  849. const uword start_slice = aux_slice1;
  850. const uword end_col_plus1 = start_col + n_cols;
  851. const uword end_row_plus1 = start_row + n_rows;
  852. const uword end_slice_plus1 = start_slice + n_slices;
  853. for(uword uslice = start_slice; uslice < end_slice_plus1; ++uslice)
  854. for(uword ucol = start_col; ucol < end_col_plus1; ++ucol )
  855. for(uword urow = start_row; urow < end_row_plus1; ++urow )
  856. {
  857. Q.at(urow, ucol, uslice) = eT( F( Q.at(urow, ucol, uslice) ) );
  858. }
  859. }
  860. //! imbue (fill) the subview with values provided by a functor
  861. template<typename eT>
  862. template<typename functor>
  863. inline
  864. void
  865. subview_cube<eT>::imbue(functor F)
  866. {
  867. arma_extra_debug_sigprint();
  868. Cube<eT>& Q = const_cast< Cube<eT>& >(m);
  869. const uword start_col = aux_col1;
  870. const uword start_row = aux_row1;
  871. const uword start_slice = aux_slice1;
  872. const uword end_col_plus1 = start_col + n_cols;
  873. const uword end_row_plus1 = start_row + n_rows;
  874. const uword end_slice_plus1 = start_slice + n_slices;
  875. for(uword uslice = start_slice; uslice < end_slice_plus1; ++uslice)
  876. for(uword ucol = start_col; ucol < end_col_plus1; ++ucol )
  877. for(uword urow = start_row; urow < end_row_plus1; ++urow )
  878. {
  879. Q.at(urow, ucol, uslice) = eT( F() );
  880. }
  881. }
  882. #if defined(ARMA_USE_CXX11)
  883. //! apply a lambda function to each slice, where each slice is interpreted as a matrix
  884. template<typename eT>
  885. inline
  886. void
  887. subview_cube<eT>::each_slice(const std::function< void(Mat<eT>&) >& F)
  888. {
  889. arma_extra_debug_sigprint();
  890. Mat<eT> tmp1(n_rows, n_cols);
  891. Mat<eT> tmp2('j', tmp1.memptr(), n_rows, n_cols);
  892. for(uword slice_id=0; slice_id < n_slices; ++slice_id)
  893. {
  894. for(uword col_id=0; col_id < n_cols; ++col_id)
  895. {
  896. arrayops::copy( tmp1.colptr(col_id), slice_colptr(slice_id, col_id), n_rows );
  897. }
  898. F(tmp2);
  899. for(uword col_id=0; col_id < n_cols; ++col_id)
  900. {
  901. arrayops::copy( slice_colptr(slice_id, col_id), tmp1.colptr(col_id), n_rows );
  902. }
  903. }
  904. }
  905. template<typename eT>
  906. inline
  907. void
  908. subview_cube<eT>::each_slice(const std::function< void(const Mat<eT>&) >& F) const
  909. {
  910. arma_extra_debug_sigprint();
  911. Mat<eT> tmp1(n_rows, n_cols);
  912. const Mat<eT> tmp2('j', tmp1.memptr(), n_rows, n_cols);
  913. for(uword slice_id=0; slice_id < n_slices; ++slice_id)
  914. {
  915. for(uword col_id=0; col_id < n_cols; ++col_id)
  916. {
  917. arrayops::copy( tmp1.colptr(col_id), slice_colptr(slice_id, col_id), n_rows );
  918. }
  919. F(tmp2);
  920. }
  921. }
  922. #endif
  923. template<typename eT>
  924. inline
  925. void
  926. subview_cube<eT>::replace(const eT old_val, const eT new_val)
  927. {
  928. arma_extra_debug_sigprint();
  929. const uword local_n_rows = n_rows;
  930. const uword local_n_cols = n_cols;
  931. const uword local_n_slices = n_slices;
  932. for(uword slice = 0; slice < local_n_slices; ++slice)
  933. {
  934. for(uword col = 0; col < local_n_cols; ++col)
  935. {
  936. arrayops::replace(slice_colptr(slice,col), local_n_rows, old_val, new_val);
  937. }
  938. }
  939. }
  940. template<typename eT>
  941. inline
  942. void
  943. subview_cube<eT>::clean(const typename get_pod_type<eT>::result threshold)
  944. {
  945. arma_extra_debug_sigprint();
  946. const uword local_n_rows = n_rows;
  947. const uword local_n_cols = n_cols;
  948. const uword local_n_slices = n_slices;
  949. for(uword slice = 0; slice < local_n_slices; ++slice)
  950. {
  951. for(uword col = 0; col < local_n_cols; ++col)
  952. {
  953. arrayops::clean( slice_colptr(slice,col), local_n_rows, threshold );
  954. }
  955. }
  956. }
  957. template<typename eT>
  958. inline
  959. void
  960. subview_cube<eT>::fill(const eT val)
  961. {
  962. arma_extra_debug_sigprint();
  963. const uword local_n_rows = n_rows;
  964. const uword local_n_cols = n_cols;
  965. const uword local_n_slices = n_slices;
  966. for(uword slice = 0; slice < local_n_slices; ++slice)
  967. {
  968. for(uword col = 0; col < local_n_cols; ++col)
  969. {
  970. arrayops::inplace_set( slice_colptr(slice,col), val, local_n_rows );
  971. }
  972. }
  973. }
  974. template<typename eT>
  975. inline
  976. void
  977. subview_cube<eT>::zeros()
  978. {
  979. arma_extra_debug_sigprint();
  980. const uword local_n_rows = n_rows;
  981. const uword local_n_cols = n_cols;
  982. const uword local_n_slices = n_slices;
  983. for(uword slice = 0; slice < local_n_slices; ++slice)
  984. {
  985. for(uword col = 0; col < local_n_cols; ++col)
  986. {
  987. arrayops::fill_zeros( slice_colptr(slice,col), local_n_rows );
  988. }
  989. }
  990. }
  991. template<typename eT>
  992. inline
  993. void
  994. subview_cube<eT>::ones()
  995. {
  996. arma_extra_debug_sigprint();
  997. fill(eT(1));
  998. }
  999. template<typename eT>
  1000. inline
  1001. void
  1002. subview_cube<eT>::randu()
  1003. {
  1004. arma_extra_debug_sigprint();
  1005. const uword local_n_rows = n_rows;
  1006. const uword local_n_cols = n_cols;
  1007. const uword local_n_slices = n_slices;
  1008. for(uword slice = 0; slice < local_n_slices; ++slice)
  1009. {
  1010. for(uword col = 0; col < local_n_cols; ++col)
  1011. {
  1012. arma_rng::randu<eT>::fill( slice_colptr(slice,col), local_n_rows );
  1013. }
  1014. }
  1015. }
  1016. template<typename eT>
  1017. inline
  1018. void
  1019. subview_cube<eT>::randn()
  1020. {
  1021. arma_extra_debug_sigprint();
  1022. const uword local_n_rows = n_rows;
  1023. const uword local_n_cols = n_cols;
  1024. const uword local_n_slices = n_slices;
  1025. for(uword slice = 0; slice < local_n_slices; ++slice)
  1026. {
  1027. for(uword col = 0; col < local_n_cols; ++col)
  1028. {
  1029. arma_rng::randn<eT>::fill( slice_colptr(slice,col), local_n_rows );
  1030. }
  1031. }
  1032. }
  1033. template<typename eT>
  1034. inline
  1035. arma_warn_unused
  1036. bool
  1037. subview_cube<eT>::is_finite() const
  1038. {
  1039. arma_extra_debug_sigprint();
  1040. const uword local_n_rows = n_rows;
  1041. const uword local_n_cols = n_cols;
  1042. const uword local_n_slices = n_slices;
  1043. for(uword slice = 0; slice < local_n_slices; ++slice)
  1044. {
  1045. for(uword col = 0; col < local_n_cols; ++col)
  1046. {
  1047. if(arrayops::is_finite(slice_colptr(slice,col), local_n_rows) == false) { return false; }
  1048. }
  1049. }
  1050. return true;
  1051. }
  1052. template<typename eT>
  1053. inline
  1054. arma_warn_unused
  1055. bool
  1056. subview_cube<eT>::is_zero(const typename get_pod_type<eT>::result tol) const
  1057. {
  1058. arma_extra_debug_sigprint();
  1059. const uword local_n_rows = n_rows;
  1060. const uword local_n_cols = n_cols;
  1061. const uword local_n_slices = n_slices;
  1062. for(uword slice = 0; slice < local_n_slices; ++slice)
  1063. {
  1064. for(uword col = 0; col < local_n_cols; ++col)
  1065. {
  1066. if(arrayops::is_zero(slice_colptr(slice,col), local_n_rows, tol) == false) { return false; }
  1067. }
  1068. }
  1069. return true;
  1070. }
  1071. template<typename eT>
  1072. inline
  1073. arma_warn_unused
  1074. bool
  1075. subview_cube<eT>::has_inf() const
  1076. {
  1077. arma_extra_debug_sigprint();
  1078. const uword local_n_rows = n_rows;
  1079. const uword local_n_cols = n_cols;
  1080. const uword local_n_slices = n_slices;
  1081. for(uword slice = 0; slice < local_n_slices; ++slice)
  1082. {
  1083. for(uword col = 0; col < local_n_cols; ++col)
  1084. {
  1085. if(arrayops::has_inf(slice_colptr(slice,col), local_n_rows)) { return true; }
  1086. }
  1087. }
  1088. return false;
  1089. }
  1090. template<typename eT>
  1091. inline
  1092. arma_warn_unused
  1093. bool
  1094. subview_cube<eT>::has_nan() const
  1095. {
  1096. arma_extra_debug_sigprint();
  1097. const uword local_n_rows = n_rows;
  1098. const uword local_n_cols = n_cols;
  1099. const uword local_n_slices = n_slices;
  1100. for(uword slice = 0; slice < local_n_slices; ++slice)
  1101. {
  1102. for(uword col = 0; col < local_n_cols; ++col)
  1103. {
  1104. if(arrayops::has_nan(slice_colptr(slice,col), local_n_rows)) { return true; }
  1105. }
  1106. }
  1107. return false;
  1108. }
  1109. template<typename eT>
  1110. inline
  1111. eT
  1112. subview_cube<eT>::at_alt(const uword i) const
  1113. {
  1114. return operator[](i);
  1115. }
  1116. template<typename eT>
  1117. inline
  1118. eT&
  1119. subview_cube<eT>::operator[](const uword i)
  1120. {
  1121. const uword in_slice = i / n_elem_slice;
  1122. const uword offset = in_slice * n_elem_slice;
  1123. const uword j = i - offset;
  1124. const uword in_col = j / n_rows;
  1125. const uword in_row = j % n_rows;
  1126. const uword index = (in_slice + aux_slice1)*m.n_elem_slice + (in_col + aux_col1)*m.n_rows + aux_row1 + in_row;
  1127. return access::rw( (const_cast< Cube<eT>& >(m)).mem[index] );
  1128. }
  1129. template<typename eT>
  1130. inline
  1131. eT
  1132. subview_cube<eT>::operator[](const uword i) const
  1133. {
  1134. const uword in_slice = i / n_elem_slice;
  1135. const uword offset = in_slice * n_elem_slice;
  1136. const uword j = i - offset;
  1137. const uword in_col = j / n_rows;
  1138. const uword in_row = j % n_rows;
  1139. const uword index = (in_slice + aux_slice1)*m.n_elem_slice + (in_col + aux_col1)*m.n_rows + aux_row1 + in_row;
  1140. return m.mem[index];
  1141. }
  1142. template<typename eT>
  1143. inline
  1144. eT&
  1145. subview_cube<eT>::operator()(const uword i)
  1146. {
  1147. arma_debug_check( (i >= n_elem), "subview_cube::operator(): index out of bounds" );
  1148. const uword in_slice = i / n_elem_slice;
  1149. const uword offset = in_slice * n_elem_slice;
  1150. const uword j = i - offset;
  1151. const uword in_col = j / n_rows;
  1152. const uword in_row = j % n_rows;
  1153. const uword index = (in_slice + aux_slice1)*m.n_elem_slice + (in_col + aux_col1)*m.n_rows + aux_row1 + in_row;
  1154. return access::rw( (const_cast< Cube<eT>& >(m)).mem[index] );
  1155. }
  1156. template<typename eT>
  1157. inline
  1158. eT
  1159. subview_cube<eT>::operator()(const uword i) const
  1160. {
  1161. arma_debug_check( (i >= n_elem), "subview_cube::operator(): index out of bounds" );
  1162. const uword in_slice = i / n_elem_slice;
  1163. const uword offset = in_slice * n_elem_slice;
  1164. const uword j = i - offset;
  1165. const uword in_col = j / n_rows;
  1166. const uword in_row = j % n_rows;
  1167. const uword index = (in_slice + aux_slice1)*m.n_elem_slice + (in_col + aux_col1)*m.n_rows + aux_row1 + in_row;
  1168. return m.mem[index];
  1169. }
  1170. template<typename eT>
  1171. arma_inline
  1172. eT&
  1173. subview_cube<eT>::operator()(const uword in_row, const uword in_col, const uword in_slice)
  1174. {
  1175. arma_debug_check( ( (in_row >= n_rows) || (in_col >= n_cols) || (in_slice >= n_slices) ), "subview_cube::operator(): location out of bounds" );
  1176. const uword index = (in_slice + aux_slice1)*m.n_elem_slice + (in_col + aux_col1)*m.n_rows + aux_row1 + in_row;
  1177. return access::rw( (const_cast< Cube<eT>& >(m)).mem[index] );
  1178. }
  1179. template<typename eT>
  1180. arma_inline
  1181. eT
  1182. subview_cube<eT>::operator()(const uword in_row, const uword in_col, const uword in_slice) const
  1183. {
  1184. arma_debug_check( ( (in_row >= n_rows) || (in_col >= n_cols) || (in_slice >= n_slices) ), "subview_cube::operator(): location out of bounds" );
  1185. const uword index = (in_slice + aux_slice1)*m.n_elem_slice + (in_col + aux_col1)*m.n_rows + aux_row1 + in_row;
  1186. return m.mem[index];
  1187. }
  1188. template<typename eT>
  1189. arma_inline
  1190. eT&
  1191. subview_cube<eT>::at(const uword in_row, const uword in_col, const uword in_slice)
  1192. {
  1193. const uword index = (in_slice + aux_slice1)*m.n_elem_slice + (in_col + aux_col1)*m.n_rows + aux_row1 + in_row;
  1194. return access::rw( (const_cast< Cube<eT>& >(m)).mem[index] );
  1195. }
  1196. template<typename eT>
  1197. arma_inline
  1198. eT
  1199. subview_cube<eT>::at(const uword in_row, const uword in_col, const uword in_slice) const
  1200. {
  1201. const uword index = (in_slice + aux_slice1)*m.n_elem_slice + (in_col + aux_col1)*m.n_rows + aux_row1 + in_row;
  1202. return m.mem[index];
  1203. }
  1204. template<typename eT>
  1205. arma_inline
  1206. eT*
  1207. subview_cube<eT>::slice_colptr(const uword in_slice, const uword in_col)
  1208. {
  1209. return & access::rw((const_cast< Cube<eT>& >(m)).mem[ (in_slice + aux_slice1)*m.n_elem_slice + (in_col + aux_col1)*m.n_rows + aux_row1 ]);
  1210. }
  1211. template<typename eT>
  1212. arma_inline
  1213. const eT*
  1214. subview_cube<eT>::slice_colptr(const uword in_slice, const uword in_col) const
  1215. {
  1216. return & m.mem[ (in_slice + aux_slice1)*m.n_elem_slice + (in_col + aux_col1)*m.n_rows + aux_row1 ];
  1217. }
  1218. template<typename eT>
  1219. inline
  1220. bool
  1221. subview_cube<eT>::check_overlap(const subview_cube<eT>& x) const
  1222. {
  1223. const subview_cube<eT>& t = *this;
  1224. if(&t.m != &x.m)
  1225. {
  1226. return false;
  1227. }
  1228. else
  1229. {
  1230. if( (t.n_elem == 0) || (x.n_elem == 0) )
  1231. {
  1232. return false;
  1233. }
  1234. else
  1235. {
  1236. const uword t_row_start = t.aux_row1;
  1237. const uword t_row_end_p1 = t_row_start + t.n_rows;
  1238. const uword t_col_start = t.aux_col1;
  1239. const uword t_col_end_p1 = t_col_start + t.n_cols;
  1240. const uword t_slice_start = t.aux_slice1;
  1241. const uword t_slice_end_p1 = t_slice_start + t.n_slices;
  1242. const uword x_row_start = x.aux_row1;
  1243. const uword x_row_end_p1 = x_row_start + x.n_rows;
  1244. const uword x_col_start = x.aux_col1;
  1245. const uword x_col_end_p1 = x_col_start + x.n_cols;
  1246. const uword x_slice_start = x.aux_slice1;
  1247. const uword x_slice_end_p1 = x_slice_start + x.n_slices;
  1248. const bool outside_rows = ( (x_row_start >= t_row_end_p1 ) || (t_row_start >= x_row_end_p1 ) );
  1249. const bool outside_cols = ( (x_col_start >= t_col_end_p1 ) || (t_col_start >= x_col_end_p1 ) );
  1250. const bool outside_slices = ( (x_slice_start >= t_slice_end_p1) || (t_slice_start >= x_slice_end_p1) );
  1251. return ( (outside_rows == false) && (outside_cols == false) && (outside_slices == false) );
  1252. }
  1253. }
  1254. }
  1255. template<typename eT>
  1256. inline
  1257. bool
  1258. subview_cube<eT>::check_overlap(const Mat<eT>& x) const
  1259. {
  1260. const subview_cube<eT>& t = *this;
  1261. const uword t_aux_slice1 = t.aux_slice1;
  1262. const uword t_aux_slice2_plus_1 = t_aux_slice1 + t.n_slices;
  1263. for(uword slice = t_aux_slice1; slice < t_aux_slice2_plus_1; ++slice)
  1264. {
  1265. if(t.m.mat_ptrs[slice] != NULL)
  1266. {
  1267. const Mat<eT>& y = *(t.m.mat_ptrs[slice]);
  1268. if( x.memptr() == y.memptr() ) { return true; }
  1269. }
  1270. }
  1271. return false;
  1272. }
  1273. //! cube X = Y.subcube(...)
  1274. template<typename eT>
  1275. inline
  1276. void
  1277. subview_cube<eT>::extract(Cube<eT>& out, const subview_cube<eT>& in)
  1278. {
  1279. arma_extra_debug_sigprint();
  1280. // NOTE: we're assuming that the cube has already been set to the correct size and there is no aliasing;
  1281. // size setting and alias checking is done by either the Cube contructor or operator=()
  1282. const uword n_rows = in.n_rows;
  1283. const uword n_cols = in.n_cols;
  1284. const uword n_slices = in.n_slices;
  1285. arma_extra_debug_print(arma_str::format("out.n_rows = %d out.n_cols = %d out.n_slices = %d in.m.n_rows = %d in.m.n_cols = %d in.m.n_slices = %d") % out.n_rows % out.n_cols % out.n_slices % in.m.n_rows % in.m.n_cols % in.m.n_slices);
  1286. for(uword slice = 0; slice < n_slices; ++slice)
  1287. {
  1288. for(uword col = 0; col < n_cols; ++col)
  1289. {
  1290. arrayops::copy( out.slice_colptr(slice,col), in.slice_colptr(slice,col), n_rows );
  1291. }
  1292. }
  1293. }
  1294. //! cube X += Y.subcube(...)
  1295. template<typename eT>
  1296. inline
  1297. void
  1298. subview_cube<eT>::plus_inplace(Cube<eT>& out, const subview_cube<eT>& in)
  1299. {
  1300. arma_extra_debug_sigprint();
  1301. arma_debug_assert_same_size(out, in, "addition");
  1302. const uword n_rows = out.n_rows;
  1303. const uword n_cols = out.n_cols;
  1304. const uword n_slices = out.n_slices;
  1305. for(uword slice = 0; slice<n_slices; ++slice)
  1306. {
  1307. for(uword col = 0; col<n_cols; ++col)
  1308. {
  1309. arrayops::inplace_plus( out.slice_colptr(slice,col), in.slice_colptr(slice,col), n_rows );
  1310. }
  1311. }
  1312. }
  1313. //! cube X -= Y.subcube(...)
  1314. template<typename eT>
  1315. inline
  1316. void
  1317. subview_cube<eT>::minus_inplace(Cube<eT>& out, const subview_cube<eT>& in)
  1318. {
  1319. arma_extra_debug_sigprint();
  1320. arma_debug_assert_same_size(out, in, "subtraction");
  1321. const uword n_rows = out.n_rows;
  1322. const uword n_cols = out.n_cols;
  1323. const uword n_slices = out.n_slices;
  1324. for(uword slice = 0; slice<n_slices; ++slice)
  1325. {
  1326. for(uword col = 0; col<n_cols; ++col)
  1327. {
  1328. arrayops::inplace_minus( out.slice_colptr(slice,col), in.slice_colptr(slice,col), n_rows );
  1329. }
  1330. }
  1331. }
  1332. //! cube X %= Y.subcube(...)
  1333. template<typename eT>
  1334. inline
  1335. void
  1336. subview_cube<eT>::schur_inplace(Cube<eT>& out, const subview_cube<eT>& in)
  1337. {
  1338. arma_extra_debug_sigprint();
  1339. arma_debug_assert_same_size(out, in, "element-wise multiplication");
  1340. const uword n_rows = out.n_rows;
  1341. const uword n_cols = out.n_cols;
  1342. const uword n_slices = out.n_slices;
  1343. for(uword slice = 0; slice<n_slices; ++slice)
  1344. {
  1345. for(uword col = 0; col<n_cols; ++col)
  1346. {
  1347. arrayops::inplace_mul( out.slice_colptr(slice,col), in.slice_colptr(slice,col), n_rows );
  1348. }
  1349. }
  1350. }
  1351. //! cube X /= Y.subcube(...)
  1352. template<typename eT>
  1353. inline
  1354. void
  1355. subview_cube<eT>::div_inplace(Cube<eT>& out, const subview_cube<eT>& in)
  1356. {
  1357. arma_extra_debug_sigprint();
  1358. arma_debug_assert_same_size(out, in, "element-wise division");
  1359. const uword n_rows = out.n_rows;
  1360. const uword n_cols = out.n_cols;
  1361. const uword n_slices = out.n_slices;
  1362. for(uword slice = 0; slice<n_slices; ++slice)
  1363. {
  1364. for(uword col = 0; col<n_cols; ++col)
  1365. {
  1366. arrayops::inplace_div( out.slice_colptr(slice,col), in.slice_colptr(slice,col), n_rows );
  1367. }
  1368. }
  1369. }
  1370. //! mat X = Y.subcube(...)
  1371. template<typename eT>
  1372. inline
  1373. void
  1374. subview_cube<eT>::extract(Mat<eT>& out, const subview_cube<eT>& in)
  1375. {
  1376. arma_extra_debug_sigprint();
  1377. arma_debug_assert_cube_as_mat(out, in, "copy into matrix", false);
  1378. const uword in_n_rows = in.n_rows;
  1379. const uword in_n_cols = in.n_cols;
  1380. const uword in_n_slices = in.n_slices;
  1381. const uword out_vec_state = out.vec_state;
  1382. if(in_n_slices == 1)
  1383. {
  1384. out.set_size(in_n_rows, in_n_cols);
  1385. for(uword col=0; col < in_n_cols; ++col)
  1386. {
  1387. arrayops::copy( out.colptr(col), in.slice_colptr(0, col), in_n_rows );
  1388. }
  1389. }
  1390. else
  1391. {
  1392. if(out_vec_state == 0)
  1393. {
  1394. if(in_n_cols == 1)
  1395. {
  1396. out.set_size(in_n_rows, in_n_slices);
  1397. for(uword i=0; i < in_n_slices; ++i)
  1398. {
  1399. arrayops::copy( out.colptr(i), in.slice_colptr(i, 0), in_n_rows );
  1400. }
  1401. }
  1402. else
  1403. if(in_n_rows == 1)
  1404. {
  1405. const Cube<eT>& Q = in.m;
  1406. const uword in_aux_row1 = in.aux_row1;
  1407. const uword in_aux_col1 = in.aux_col1;
  1408. const uword in_aux_slice1 = in.aux_slice1;
  1409. out.set_size(in_n_cols, in_n_slices);
  1410. for(uword slice=0; slice < in_n_slices; ++slice)
  1411. {
  1412. const uword mod_slice = in_aux_slice1 + slice;
  1413. eT* out_colptr = out.colptr(slice);
  1414. uword i,j;
  1415. for(i=0, j=1; j < in_n_cols; i+=2, j+=2)
  1416. {
  1417. const eT tmp_i = Q.at(in_aux_row1, in_aux_col1 + i, mod_slice);
  1418. const eT tmp_j = Q.at(in_aux_row1, in_aux_col1 + j, mod_slice);
  1419. out_colptr[i] = tmp_i;
  1420. out_colptr[j] = tmp_j;
  1421. }
  1422. if(i < in_n_cols)
  1423. {
  1424. out_colptr[i] = Q.at(in_aux_row1, in_aux_col1 + i, mod_slice);
  1425. }
  1426. }
  1427. }
  1428. }
  1429. else
  1430. {
  1431. out.set_size(in_n_slices);
  1432. eT* out_mem = out.memptr();
  1433. const Cube<eT>& Q = in.m;
  1434. const uword in_aux_row1 = in.aux_row1;
  1435. const uword in_aux_col1 = in.aux_col1;
  1436. const uword in_aux_slice1 = in.aux_slice1;
  1437. for(uword i=0; i<in_n_slices; ++i)
  1438. {
  1439. out_mem[i] = Q.at(in_aux_row1, in_aux_col1, in_aux_slice1 + i);
  1440. }
  1441. }
  1442. }
  1443. }
  1444. //! mat X += Y.subcube(...)
  1445. template<typename eT>
  1446. inline
  1447. void
  1448. subview_cube<eT>::plus_inplace(Mat<eT>& out, const subview_cube<eT>& in)
  1449. {
  1450. arma_extra_debug_sigprint();
  1451. arma_debug_assert_cube_as_mat(out, in, "addition", true);
  1452. const uword in_n_rows = in.n_rows;
  1453. const uword in_n_cols = in.n_cols;
  1454. const uword in_n_slices = in.n_slices;
  1455. const uword out_n_rows = out.n_rows;
  1456. const uword out_n_cols = out.n_cols;
  1457. const uword out_vec_state = out.vec_state;
  1458. if(in_n_slices == 1)
  1459. {
  1460. if( (arma_config::debug) && ((out_n_rows != in_n_rows) || (out_n_cols != in_n_cols)) )
  1461. {
  1462. std::ostringstream tmp;
  1463. tmp
  1464. << "in-place addition: "
  1465. << out_n_rows << 'x' << out_n_cols << " output matrix is incompatible with "
  1466. << in_n_rows << 'x' << in_n_cols << 'x' << in_n_slices << " cube interpreted as "
  1467. << in_n_rows << 'x' << in_n_cols << " matrix";
  1468. arma_stop_logic_error(tmp.str());
  1469. }
  1470. for(uword col=0; col < in_n_cols; ++col)
  1471. {
  1472. arrayops::inplace_plus( out.colptr(col), in.slice_colptr(0, col), in_n_rows );
  1473. }
  1474. }
  1475. else
  1476. {
  1477. if(out_vec_state == 0)
  1478. {
  1479. if( (in_n_rows == out_n_rows) && (in_n_cols == 1) && (in_n_slices == out_n_cols) )
  1480. {
  1481. for(uword i=0; i < in_n_slices; ++i)
  1482. {
  1483. arrayops::inplace_plus( out.colptr(i), in.slice_colptr(i, 0), in_n_rows );
  1484. }
  1485. }
  1486. else
  1487. if( (in_n_rows == 1) && (in_n_cols == out_n_rows) && (in_n_slices == out_n_cols) )
  1488. {
  1489. const Cube<eT>& Q = in.m;
  1490. const uword in_aux_row1 = in.aux_row1;
  1491. const uword in_aux_col1 = in.aux_col1;
  1492. const uword in_aux_slice1 = in.aux_slice1;
  1493. for(uword slice=0; slice < in_n_slices; ++slice)
  1494. {
  1495. const uword mod_slice = in_aux_slice1 + slice;
  1496. eT* out_colptr = out.colptr(slice);
  1497. uword i,j;
  1498. for(i=0, j=1; j < in_n_cols; i+=2, j+=2)
  1499. {
  1500. const eT tmp_i = Q.at(in_aux_row1, in_aux_col1 + i, mod_slice);
  1501. const eT tmp_j = Q.at(in_aux_row1, in_aux_col1 + j, mod_slice);
  1502. out_colptr[i] += tmp_i;
  1503. out_colptr[j] += tmp_j;
  1504. }
  1505. if(i < in_n_cols)
  1506. {
  1507. out_colptr[i] += Q.at(in_aux_row1, in_aux_col1 + i, mod_slice);
  1508. }
  1509. }
  1510. }
  1511. }
  1512. else
  1513. {
  1514. eT* out_mem = out.memptr();
  1515. const Cube<eT>& Q = in.m;
  1516. const uword in_aux_row1 = in.aux_row1;
  1517. const uword in_aux_col1 = in.aux_col1;
  1518. const uword in_aux_slice1 = in.aux_slice1;
  1519. for(uword i=0; i<in_n_slices; ++i)
  1520. {
  1521. out_mem[i] += Q.at(in_aux_row1, in_aux_col1, in_aux_slice1 + i);
  1522. }
  1523. }
  1524. }
  1525. }
  1526. //! mat X -= Y.subcube(...)
  1527. template<typename eT>
  1528. inline
  1529. void
  1530. subview_cube<eT>::minus_inplace(Mat<eT>& out, const subview_cube<eT>& in)
  1531. {
  1532. arma_extra_debug_sigprint();
  1533. arma_debug_assert_cube_as_mat(out, in, "subtraction", true);
  1534. const uword in_n_rows = in.n_rows;
  1535. const uword in_n_cols = in.n_cols;
  1536. const uword in_n_slices = in.n_slices;
  1537. const uword out_n_rows = out.n_rows;
  1538. const uword out_n_cols = out.n_cols;
  1539. const uword out_vec_state = out.vec_state;
  1540. if(in_n_slices == 1)
  1541. {
  1542. if( (arma_config::debug) && ((out_n_rows != in_n_rows) || (out_n_cols != in_n_cols)) )
  1543. {
  1544. std::ostringstream tmp;
  1545. tmp
  1546. << "in-place subtraction: "
  1547. << out_n_rows << 'x' << out_n_cols << " output matrix is incompatible with "
  1548. << in_n_rows << 'x' << in_n_cols << 'x' << in_n_slices << " cube interpreted as "
  1549. << in_n_rows << 'x' << in_n_cols << " matrix";
  1550. arma_stop_logic_error(tmp.str());
  1551. }
  1552. for(uword col=0; col < in_n_cols; ++col)
  1553. {
  1554. arrayops::inplace_minus( out.colptr(col), in.slice_colptr(0, col), in_n_rows );
  1555. }
  1556. }
  1557. else
  1558. {
  1559. if(out_vec_state == 0)
  1560. {
  1561. if( (in_n_rows == out_n_rows) && (in_n_cols == 1) && (in_n_slices == out_n_cols) )
  1562. {
  1563. for(uword i=0; i < in_n_slices; ++i)
  1564. {
  1565. arrayops::inplace_minus( out.colptr(i), in.slice_colptr(i, 0), in_n_rows );
  1566. }
  1567. }
  1568. else
  1569. if( (in_n_rows == 1) && (in_n_cols == out_n_rows) && (in_n_slices == out_n_cols) )
  1570. {
  1571. const Cube<eT>& Q = in.m;
  1572. const uword in_aux_row1 = in.aux_row1;
  1573. const uword in_aux_col1 = in.aux_col1;
  1574. const uword in_aux_slice1 = in.aux_slice1;
  1575. for(uword slice=0; slice < in_n_slices; ++slice)
  1576. {
  1577. const uword mod_slice = in_aux_slice1 + slice;
  1578. eT* out_colptr = out.colptr(slice);
  1579. uword i,j;
  1580. for(i=0, j=1; j < in_n_cols; i+=2, j+=2)
  1581. {
  1582. const eT tmp_i = Q.at(in_aux_row1, in_aux_col1 + i, mod_slice);
  1583. const eT tmp_j = Q.at(in_aux_row1, in_aux_col1 + j, mod_slice);
  1584. out_colptr[i] -= tmp_i;
  1585. out_colptr[j] -= tmp_j;
  1586. }
  1587. if(i < in_n_cols)
  1588. {
  1589. out_colptr[i] -= Q.at(in_aux_row1, in_aux_col1 + i, mod_slice);
  1590. }
  1591. }
  1592. }
  1593. }
  1594. else
  1595. {
  1596. eT* out_mem = out.memptr();
  1597. const Cube<eT>& Q = in.m;
  1598. const uword in_aux_row1 = in.aux_row1;
  1599. const uword in_aux_col1 = in.aux_col1;
  1600. const uword in_aux_slice1 = in.aux_slice1;
  1601. for(uword i=0; i<in_n_slices; ++i)
  1602. {
  1603. out_mem[i] -= Q.at(in_aux_row1, in_aux_col1, in_aux_slice1 + i);
  1604. }
  1605. }
  1606. }
  1607. }
  1608. //! mat X %= Y.subcube(...)
  1609. template<typename eT>
  1610. inline
  1611. void
  1612. subview_cube<eT>::schur_inplace(Mat<eT>& out, const subview_cube<eT>& in)
  1613. {
  1614. arma_extra_debug_sigprint();
  1615. arma_debug_assert_cube_as_mat(out, in, "element-wise multiplication", true);
  1616. const uword in_n_rows = in.n_rows;
  1617. const uword in_n_cols = in.n_cols;
  1618. const uword in_n_slices = in.n_slices;
  1619. const uword out_n_rows = out.n_rows;
  1620. const uword out_n_cols = out.n_cols;
  1621. const uword out_vec_state = out.vec_state;
  1622. if(in_n_slices == 1)
  1623. {
  1624. if( (arma_config::debug) && ((out_n_rows != in_n_rows) || (out_n_cols != in_n_cols)) )
  1625. {
  1626. std::ostringstream tmp;
  1627. tmp
  1628. << "in-place element-wise multiplication: "
  1629. << out_n_rows << 'x' << out_n_cols << " output matrix is incompatible with "
  1630. << in_n_rows << 'x' << in_n_cols << 'x' << in_n_slices << " cube interpreted as "
  1631. << in_n_rows << 'x' << in_n_cols << " matrix";
  1632. arma_stop_logic_error(tmp.str());
  1633. }
  1634. for(uword col=0; col < in_n_cols; ++col)
  1635. {
  1636. arrayops::inplace_mul( out.colptr(col), in.slice_colptr(0, col), in_n_rows );
  1637. }
  1638. }
  1639. else
  1640. {
  1641. if(out_vec_state == 0)
  1642. {
  1643. if( (in_n_rows == out_n_rows) && (in_n_cols == 1) && (in_n_slices == out_n_cols) )
  1644. {
  1645. for(uword i=0; i < in_n_slices; ++i)
  1646. {
  1647. arrayops::inplace_mul( out.colptr(i), in.slice_colptr(i, 0), in_n_rows );
  1648. }
  1649. }
  1650. else
  1651. if( (in_n_rows == 1) && (in_n_cols == out_n_rows) && (in_n_slices == out_n_cols) )
  1652. {
  1653. const Cube<eT>& Q = in.m;
  1654. const uword in_aux_row1 = in.aux_row1;
  1655. const uword in_aux_col1 = in.aux_col1;
  1656. const uword in_aux_slice1 = in.aux_slice1;
  1657. for(uword slice=0; slice < in_n_slices; ++slice)
  1658. {
  1659. const uword mod_slice = in_aux_slice1 + slice;
  1660. eT* out_colptr = out.colptr(slice);
  1661. uword i,j;
  1662. for(i=0, j=1; j < in_n_cols; i+=2, j+=2)
  1663. {
  1664. const eT tmp_i = Q.at(in_aux_row1, in_aux_col1 + i, mod_slice);
  1665. const eT tmp_j = Q.at(in_aux_row1, in_aux_col1 + j, mod_slice);
  1666. out_colptr[i] *= tmp_i;
  1667. out_colptr[j] *= tmp_j;
  1668. }
  1669. if(i < in_n_cols)
  1670. {
  1671. out_colptr[i] *= Q.at(in_aux_row1, in_aux_col1 + i, mod_slice);
  1672. }
  1673. }
  1674. }
  1675. }
  1676. else
  1677. {
  1678. eT* out_mem = out.memptr();
  1679. const Cube<eT>& Q = in.m;
  1680. const uword in_aux_row1 = in.aux_row1;
  1681. const uword in_aux_col1 = in.aux_col1;
  1682. const uword in_aux_slice1 = in.aux_slice1;
  1683. for(uword i=0; i<in_n_slices; ++i)
  1684. {
  1685. out_mem[i] *= Q.at(in_aux_row1, in_aux_col1, in_aux_slice1 + i);
  1686. }
  1687. }
  1688. }
  1689. }
  1690. //! mat X /= Y.subcube(...)
  1691. template<typename eT>
  1692. inline
  1693. void
  1694. subview_cube<eT>::div_inplace(Mat<eT>& out, const subview_cube<eT>& in)
  1695. {
  1696. arma_extra_debug_sigprint();
  1697. arma_debug_assert_cube_as_mat(out, in, "element-wise division", true);
  1698. const uword in_n_rows = in.n_rows;
  1699. const uword in_n_cols = in.n_cols;
  1700. const uword in_n_slices = in.n_slices;
  1701. const uword out_n_rows = out.n_rows;
  1702. const uword out_n_cols = out.n_cols;
  1703. const uword out_vec_state = out.vec_state;
  1704. if(in_n_slices == 1)
  1705. {
  1706. if( (arma_config::debug) && ((out_n_rows != in_n_rows) || (out_n_cols != in_n_cols)) )
  1707. {
  1708. std::ostringstream tmp;
  1709. tmp
  1710. << "in-place element-wise division: "
  1711. << out_n_rows << 'x' << out_n_cols << " output matrix is incompatible with "
  1712. << in_n_rows << 'x' << in_n_cols << 'x' << in_n_slices << " cube interpreted as "
  1713. << in_n_rows << 'x' << in_n_cols << " matrix";
  1714. arma_stop_logic_error(tmp.str());
  1715. }
  1716. for(uword col=0; col < in_n_cols; ++col)
  1717. {
  1718. arrayops::inplace_div( out.colptr(col), in.slice_colptr(0, col), in_n_rows );
  1719. }
  1720. }
  1721. else
  1722. {
  1723. if(out_vec_state == 0)
  1724. {
  1725. if( (in_n_rows == out_n_rows) && (in_n_cols == 1) && (in_n_slices == out_n_cols) )
  1726. {
  1727. for(uword i=0; i < in_n_slices; ++i)
  1728. {
  1729. arrayops::inplace_div( out.colptr(i), in.slice_colptr(i, 0), in_n_rows );
  1730. }
  1731. }
  1732. else
  1733. if( (in_n_rows == 1) && (in_n_cols == out_n_rows) && (in_n_slices == out_n_cols) )
  1734. {
  1735. const Cube<eT>& Q = in.m;
  1736. const uword in_aux_row1 = in.aux_row1;
  1737. const uword in_aux_col1 = in.aux_col1;
  1738. const uword in_aux_slice1 = in.aux_slice1;
  1739. for(uword slice=0; slice < in_n_slices; ++slice)
  1740. {
  1741. const uword mod_slice = in_aux_slice1 + slice;
  1742. eT* out_colptr = out.colptr(slice);
  1743. uword i,j;
  1744. for(i=0, j=1; j < in_n_cols; i+=2, j+=2)
  1745. {
  1746. const eT tmp_i = Q.at(in_aux_row1, in_aux_col1 + i, mod_slice);
  1747. const eT tmp_j = Q.at(in_aux_row1, in_aux_col1 + j, mod_slice);
  1748. out_colptr[i] /= tmp_i;
  1749. out_colptr[j] /= tmp_j;
  1750. }
  1751. if(i < in_n_cols)
  1752. {
  1753. out_colptr[i] /= Q.at(in_aux_row1, in_aux_col1 + i, mod_slice);
  1754. }
  1755. }
  1756. }
  1757. }
  1758. else
  1759. {
  1760. eT* out_mem = out.memptr();
  1761. const Cube<eT>& Q = in.m;
  1762. const uword in_aux_row1 = in.aux_row1;
  1763. const uword in_aux_col1 = in.aux_col1;
  1764. const uword in_aux_slice1 = in.aux_slice1;
  1765. for(uword i=0; i<in_n_slices; ++i)
  1766. {
  1767. out_mem[i] /= Q.at(in_aux_row1, in_aux_col1, in_aux_slice1 + i);
  1768. }
  1769. }
  1770. }
  1771. }
  1772. template<typename eT>
  1773. inline
  1774. typename subview_cube<eT>::iterator
  1775. subview_cube<eT>::begin()
  1776. {
  1777. return iterator(*this, aux_row1, aux_col1, aux_slice1);
  1778. }
  1779. template<typename eT>
  1780. inline
  1781. typename subview_cube<eT>::const_iterator
  1782. subview_cube<eT>::begin() const
  1783. {
  1784. return const_iterator(*this, aux_row1, aux_col1, aux_slice1);
  1785. }
  1786. template<typename eT>
  1787. inline
  1788. typename subview_cube<eT>::const_iterator
  1789. subview_cube<eT>::cbegin() const
  1790. {
  1791. return const_iterator(*this, aux_row1, aux_col1, aux_slice1);
  1792. }
  1793. template<typename eT>
  1794. inline
  1795. typename subview_cube<eT>::iterator
  1796. subview_cube<eT>::end()
  1797. {
  1798. return iterator(*this, aux_row1, aux_col1, aux_slice1 + n_slices);
  1799. }
  1800. template<typename eT>
  1801. inline
  1802. typename subview_cube<eT>::const_iterator
  1803. subview_cube<eT>::end() const
  1804. {
  1805. return const_iterator(*this, aux_row1, aux_col1, aux_slice1 + n_slices);
  1806. }
  1807. template<typename eT>
  1808. inline
  1809. typename subview_cube<eT>::const_iterator
  1810. subview_cube<eT>::cend() const
  1811. {
  1812. return const_iterator(*this, aux_row1, aux_col1, aux_slice1 + n_slices);
  1813. }
  1814. //
  1815. //
  1816. //
  1817. template<typename eT>
  1818. inline
  1819. subview_cube<eT>::iterator::iterator()
  1820. : M (NULL)
  1821. , current_ptr (NULL)
  1822. , current_row (0 )
  1823. , current_col (0 )
  1824. , current_slice(0 )
  1825. , aux_row1 (0 )
  1826. , aux_col1 (0 )
  1827. , aux_row2_p1 (0 )
  1828. , aux_col2_p1 (0 )
  1829. {
  1830. arma_extra_debug_sigprint();
  1831. // Technically this iterator is invalid (it does not point to a valid element)
  1832. }
  1833. template<typename eT>
  1834. inline
  1835. subview_cube<eT>::iterator::iterator(const iterator& X)
  1836. : M (X.M )
  1837. , current_ptr (X.current_ptr )
  1838. , current_row (X.current_row )
  1839. , current_col (X.current_col )
  1840. , current_slice(X.current_slice)
  1841. , aux_row1 (X.aux_row1 )
  1842. , aux_col1 (X.aux_col1 )
  1843. , aux_row2_p1 (X.aux_row2_p1 )
  1844. , aux_col2_p1 (X.aux_col2_p1 )
  1845. {
  1846. arma_extra_debug_sigprint();
  1847. }
  1848. template<typename eT>
  1849. inline
  1850. subview_cube<eT>::iterator::iterator(subview_cube<eT>& in_sv, const uword in_row, const uword in_col, const uword in_slice)
  1851. : M (&(const_cast< Cube<eT>& >(in_sv.m)))
  1852. , current_ptr (&(M->at(in_row,in_col,in_slice)) )
  1853. , current_row (in_row )
  1854. , current_col (in_col )
  1855. , current_slice(in_slice )
  1856. , aux_row1 (in_sv.aux_row1 )
  1857. , aux_col1 (in_sv.aux_col1 )
  1858. , aux_row2_p1 (in_sv.aux_row1 + in_sv.n_rows )
  1859. , aux_col2_p1 (in_sv.aux_col1 + in_sv.n_cols )
  1860. {
  1861. arma_extra_debug_sigprint();
  1862. }
  1863. template<typename eT>
  1864. inline
  1865. arma_warn_unused
  1866. eT&
  1867. subview_cube<eT>::iterator::operator*()
  1868. {
  1869. return (*current_ptr);
  1870. }
  1871. template<typename eT>
  1872. inline
  1873. typename subview_cube<eT>::iterator&
  1874. subview_cube<eT>::iterator::operator++()
  1875. {
  1876. current_row++;
  1877. if(current_row == aux_row2_p1)
  1878. {
  1879. current_row = aux_row1;
  1880. current_col++;
  1881. if(current_col == aux_col2_p1)
  1882. {
  1883. current_col = aux_col1;
  1884. current_slice++;
  1885. }
  1886. current_ptr = &( (*M).at(current_row,current_col,current_slice) );
  1887. }
  1888. else
  1889. {
  1890. current_ptr++;
  1891. }
  1892. return *this;
  1893. }
  1894. template<typename eT>
  1895. inline
  1896. arma_warn_unused
  1897. typename subview_cube<eT>::iterator
  1898. subview_cube<eT>::iterator::operator++(int)
  1899. {
  1900. typename subview_cube<eT>::iterator temp(*this);
  1901. ++(*this);
  1902. return temp;
  1903. }
  1904. template<typename eT>
  1905. inline
  1906. arma_warn_unused
  1907. bool
  1908. subview_cube<eT>::iterator::operator==(const iterator& rhs) const
  1909. {
  1910. return (current_ptr == rhs.current_ptr);
  1911. }
  1912. template<typename eT>
  1913. inline
  1914. arma_warn_unused
  1915. bool
  1916. subview_cube<eT>::iterator::operator!=(const iterator& rhs) const
  1917. {
  1918. return (current_ptr != rhs.current_ptr);
  1919. }
  1920. template<typename eT>
  1921. inline
  1922. arma_warn_unused
  1923. bool
  1924. subview_cube<eT>::iterator::operator==(const const_iterator& rhs) const
  1925. {
  1926. return (current_ptr == rhs.current_ptr);
  1927. }
  1928. template<typename eT>
  1929. inline
  1930. arma_warn_unused
  1931. bool
  1932. subview_cube<eT>::iterator::operator!=(const const_iterator& rhs) const
  1933. {
  1934. return (current_ptr != rhs.current_ptr);
  1935. }
  1936. //
  1937. //
  1938. //
  1939. template<typename eT>
  1940. inline
  1941. subview_cube<eT>::const_iterator::const_iterator()
  1942. : M (NULL)
  1943. , current_ptr (NULL)
  1944. , current_row (0 )
  1945. , current_col (0 )
  1946. , current_slice(0 )
  1947. , aux_row1 (0 )
  1948. , aux_col1 (0 )
  1949. , aux_row2_p1 (0 )
  1950. , aux_col2_p1 (0 )
  1951. {
  1952. arma_extra_debug_sigprint();
  1953. // Technically this iterator is invalid (it does not point to a valid element)
  1954. }
  1955. template<typename eT>
  1956. inline
  1957. subview_cube<eT>::const_iterator::const_iterator(const iterator& X)
  1958. : M (X.M )
  1959. , current_ptr (X.current_ptr )
  1960. , current_row (X.current_row )
  1961. , current_col (X.current_col )
  1962. , current_slice(X.current_slice)
  1963. , aux_row1 (X.aux_row1 )
  1964. , aux_col1 (X.aux_col1 )
  1965. , aux_row2_p1 (X.aux_row2_p1 )
  1966. , aux_col2_p1 (X.aux_col2_p1 )
  1967. {
  1968. arma_extra_debug_sigprint();
  1969. }
  1970. template<typename eT>
  1971. inline
  1972. subview_cube<eT>::const_iterator::const_iterator(const const_iterator& X)
  1973. : M (X.M )
  1974. , current_ptr (X.current_ptr )
  1975. , current_row (X.current_row )
  1976. , current_col (X.current_col )
  1977. , current_slice(X.current_slice)
  1978. , aux_row1 (X.aux_row1 )
  1979. , aux_col1 (X.aux_col1 )
  1980. , aux_row2_p1 (X.aux_row2_p1 )
  1981. , aux_col2_p1 (X.aux_col2_p1 )
  1982. {
  1983. arma_extra_debug_sigprint();
  1984. }
  1985. template<typename eT>
  1986. inline
  1987. subview_cube<eT>::const_iterator::const_iterator(const subview_cube<eT>& in_sv, const uword in_row, const uword in_col, const uword in_slice)
  1988. : M (&(in_sv.m) )
  1989. , current_ptr (&(M->at(in_row,in_col,in_slice)))
  1990. , current_row (in_row )
  1991. , current_col (in_col )
  1992. , current_slice(in_slice )
  1993. , aux_row1 (in_sv.aux_row1 )
  1994. , aux_col1 (in_sv.aux_col1 )
  1995. , aux_row2_p1 (in_sv.aux_row1 + in_sv.n_rows )
  1996. , aux_col2_p1 (in_sv.aux_col1 + in_sv.n_cols )
  1997. {
  1998. arma_extra_debug_sigprint();
  1999. }
  2000. template<typename eT>
  2001. inline
  2002. arma_warn_unused
  2003. const eT&
  2004. subview_cube<eT>::const_iterator::operator*()
  2005. {
  2006. return (*current_ptr);
  2007. }
  2008. template<typename eT>
  2009. inline
  2010. typename subview_cube<eT>::const_iterator&
  2011. subview_cube<eT>::const_iterator::operator++()
  2012. {
  2013. current_row++;
  2014. if(current_row == aux_row2_p1)
  2015. {
  2016. current_row = aux_row1;
  2017. current_col++;
  2018. if(current_col == aux_col2_p1)
  2019. {
  2020. current_col = aux_col1;
  2021. current_slice++;
  2022. }
  2023. current_ptr = &( (*M).at(current_row,current_col,current_slice) );
  2024. }
  2025. else
  2026. {
  2027. current_ptr++;
  2028. }
  2029. return *this;
  2030. }
  2031. template<typename eT>
  2032. inline
  2033. arma_warn_unused
  2034. typename subview_cube<eT>::const_iterator
  2035. subview_cube<eT>::const_iterator::operator++(int)
  2036. {
  2037. typename subview_cube<eT>::const_iterator temp(*this);
  2038. ++(*this);
  2039. return temp;
  2040. }
  2041. template<typename eT>
  2042. inline
  2043. arma_warn_unused
  2044. bool
  2045. subview_cube<eT>::const_iterator::operator==(const iterator& rhs) const
  2046. {
  2047. return (current_ptr == rhs.current_ptr);
  2048. }
  2049. template<typename eT>
  2050. inline
  2051. arma_warn_unused
  2052. bool
  2053. subview_cube<eT>::const_iterator::operator!=(const iterator& rhs) const
  2054. {
  2055. return (current_ptr != rhs.current_ptr);
  2056. }
  2057. template<typename eT>
  2058. inline
  2059. arma_warn_unused
  2060. bool
  2061. subview_cube<eT>::const_iterator::operator==(const const_iterator& rhs) const
  2062. {
  2063. return (current_ptr == rhs.current_ptr);
  2064. }
  2065. template<typename eT>
  2066. inline
  2067. arma_warn_unused
  2068. bool
  2069. subview_cube<eT>::const_iterator::operator!=(const const_iterator& rhs) const
  2070. {
  2071. return (current_ptr != rhs.current_ptr);
  2072. }
  2073. //! @}