mkl_direct_lapack.h 113 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920921922923924925926927928929930931932933934935936937938939940941942943944945946947948949950951952953954955956957958959960961962963964965966967968969970971972973974975976977978979980981982983984985986987988989990991992993994995996997998999100010011002100310041005100610071008100910101011101210131014101510161017101810191020102110221023102410251026102710281029103010311032103310341035103610371038103910401041104210431044104510461047104810491050105110521053105410551056105710581059106010611062106310641065106610671068106910701071107210731074107510761077107810791080108110821083108410851086108710881089109010911092109310941095109610971098109911001101110211031104110511061107110811091110111111121113111411151116111711181119112011211122112311241125112611271128112911301131113211331134113511361137113811391140114111421143114411451146114711481149115011511152115311541155115611571158115911601161116211631164116511661167116811691170117111721173117411751176117711781179118011811182118311841185118611871188118911901191119211931194119511961197119811991200120112021203120412051206120712081209121012111212121312141215121612171218121912201221122212231224122512261227122812291230123112321233123412351236123712381239124012411242124312441245124612471248124912501251125212531254125512561257125812591260126112621263126412651266126712681269127012711272127312741275127612771278127912801281128212831284128512861287128812891290129112921293129412951296129712981299130013011302130313041305130613071308130913101311131213131314131513161317131813191320132113221323132413251326132713281329133013311332133313341335133613371338133913401341134213431344134513461347134813491350135113521353135413551356135713581359136013611362136313641365136613671368136913701371137213731374137513761377137813791380138113821383138413851386138713881389139013911392139313941395139613971398139914001401140214031404140514061407140814091410141114121413141414151416141714181419142014211422142314241425142614271428142914301431143214331434143514361437143814391440144114421443144414451446144714481449145014511452145314541455145614571458145914601461146214631464146514661467146814691470147114721473147414751476147714781479148014811482148314841485148614871488148914901491149214931494149514961497149814991500150115021503150415051506150715081509151015111512151315141515151615171518151915201521152215231524152515261527152815291530153115321533153415351536153715381539154015411542154315441545154615471548154915501551155215531554155515561557155815591560156115621563156415651566156715681569157015711572157315741575157615771578157915801581158215831584158515861587158815891590159115921593159415951596159715981599160016011602160316041605160616071608160916101611161216131614161516161617161816191620162116221623162416251626162716281629163016311632163316341635163616371638163916401641164216431644164516461647164816491650165116521653165416551656165716581659166016611662166316641665166616671668166916701671167216731674167516761677167816791680168116821683168416851686168716881689169016911692169316941695169616971698169917001701170217031704170517061707170817091710171117121713171417151716171717181719172017211722172317241725172617271728172917301731173217331734173517361737173817391740174117421743174417451746174717481749175017511752175317541755175617571758175917601761176217631764176517661767176817691770177117721773177417751776177717781779178017811782178317841785178617871788178917901791179217931794179517961797179817991800180118021803180418051806180718081809181018111812181318141815181618171818181918201821182218231824182518261827182818291830183118321833183418351836183718381839184018411842184318441845184618471848184918501851185218531854185518561857185818591860186118621863186418651866186718681869187018711872187318741875187618771878187918801881188218831884188518861887188818891890189118921893189418951896189718981899190019011902190319041905190619071908190919101911191219131914191519161917191819191920192119221923192419251926192719281929193019311932193319341935193619371938193919401941194219431944194519461947194819491950195119521953195419551956195719581959196019611962196319641965196619671968196919701971197219731974197519761977197819791980198119821983198419851986198719881989199019911992199319941995199619971998199920002001200220032004200520062007200820092010201120122013201420152016201720182019202020212022202320242025202620272028202920302031203220332034203520362037203820392040204120422043204420452046204720482049205020512052205320542055205620572058205920602061206220632064206520662067206820692070207120722073207420752076207720782079208020812082208320842085208620872088208920902091209220932094209520962097209820992100210121022103210421052106210721082109211021112112211321142115211621172118211921202121212221232124212521262127212821292130213121322133213421352136213721382139214021412142214321442145214621472148214921502151215221532154215521562157215821592160216121622163216421652166216721682169217021712172217321742175217621772178217921802181218221832184218521862187218821892190219121922193219421952196219721982199220022012202220322042205220622072208220922102211221222132214221522162217221822192220222122222223222422252226222722282229223022312232223322342235223622372238223922402241224222432244224522462247224822492250225122522253225422552256225722582259226022612262226322642265226622672268226922702271227222732274227522762277227822792280228122822283228422852286228722882289229022912292229322942295229622972298229923002301230223032304230523062307230823092310231123122313231423152316231723182319232023212322232323242325232623272328232923302331233223332334233523362337233823392340234123422343234423452346234723482349235023512352235323542355235623572358235923602361236223632364236523662367236823692370237123722373237423752376237723782379238023812382238323842385238623872388238923902391239223932394239523962397239823992400240124022403240424052406240724082409241024112412241324142415241624172418241924202421242224232424242524262427242824292430243124322433243424352436243724382439244024412442244324442445244624472448244924502451245224532454245524562457245824592460246124622463246424652466246724682469247024712472247324742475247624772478247924802481248224832484248524862487248824892490249124922493249424952496249724982499250025012502250325042505250625072508250925102511251225132514251525162517251825192520252125222523252425252526252725282529253025312532253325342535253625372538253925402541254225432544254525462547254825492550255125522553255425552556255725582559256025612562256325642565256625672568256925702571257225732574257525762577257825792580258125822583258425852586258725882589259025912592259325942595259625972598259926002601260226032604260526062607260826092610261126122613261426152616261726182619262026212622262326242625262626272628262926302631263226332634263526362637263826392640264126422643264426452646264726482649265026512652265326542655265626572658265926602661266226632664266526662667266826692670267126722673267426752676267726782679268026812682268326842685268626872688268926902691269226932694269526962697269826992700270127022703270427052706270727082709271027112712271327142715271627172718271927202721272227232724272527262727272827292730273127322733273427352736273727382739274027412742274327442745274627472748274927502751275227532754275527562757275827592760276127622763276427652766276727682769277027712772277327742775277627772778277927802781278227832784278527862787278827892790279127922793279427952796279727982799280028012802280328042805280628072808280928102811281228132814281528162817281828192820282128222823282428252826282728282829283028312832283328342835283628372838283928402841284228432844284528462847284828492850285128522853285428552856285728582859286028612862286328642865286628672868286928702871287228732874287528762877287828792880288128822883288428852886288728882889289028912892289328942895289628972898289929002901290229032904290529062907290829092910291129122913291429152916291729182919292029212922292329242925292629272928292929302931293229332934293529362937293829392940294129422943294429452946294729482949295029512952295329542955295629572958295929602961296229632964296529662967296829692970297129722973297429752976297729782979298029812982298329842985298629872988298929902991299229932994299529962997299829993000300130023003300430053006300730083009301030113012301330143015301630173018301930203021302230233024302530263027302830293030303130323033303430353036303730383039304030413042304330443045304630473048304930503051305230533054305530563057305830593060306130623063306430653066306730683069307030713072307330743075307630773078307930803081308230833084308530863087308830893090309130923093309430953096309730983099310031013102310331043105310631073108310931103111311231133114311531163117311831193120312131223123312431253126312731283129313031313132313331343135313631373138313931403141314231433144314531463147314831493150315131523153315431553156315731583159316031613162316331643165316631673168316931703171317231733174317531763177317831793180318131823183318431853186318731883189319031913192319331943195319631973198319932003201320232033204320532063207320832093210321132123213321432153216321732183219322032213222322332243225322632273228322932303231323232333234323532363237323832393240324132423243324432453246324732483249325032513252325332543255325632573258325932603261
  1. /*******************************************************************************
  2. * Copyright 2016-2022 Intel Corporation.
  3. *
  4. * This software and the related documents are Intel copyrighted materials, and
  5. * your use of them is governed by the express license under which they were
  6. * provided to you (License). Unless the License provides otherwise, you may not
  7. * use, modify, copy, publish, distribute, disclose or transmit this software or
  8. * the related documents without Intel's prior written permission.
  9. *
  10. * This software and the related documents are provided as is, with no express
  11. * or implied warranties, other than those that are expressly stated in the
  12. * License.
  13. *******************************************************************************/
  14. /*
  15. ! Content:
  16. ! Intel(R) oneAPI Math Kernel Library (oneMKL) C functions that can be inlined
  17. !******************************************************************************/
  18. #ifndef MKL_DC_UNSAFE
  19. #define MKL_DC_UNSAFE 0
  20. #endif
  21. #undef mkl_dc_getrf
  22. #undef mkl_dc_lapacke_getrf_convert
  23. #undef mkl_dc_getrf_1x1
  24. #undef mkl_dc_getrf_2x2
  25. #undef mkl_dc_getrf_3x3
  26. #undef mkl_dc_getrf_4x4
  27. #undef mkl_dc_getrf_5x5
  28. #undef mkl_dc_getrf_mxn
  29. #undef mkl_dc_getrf_generic
  30. #undef mkl_dc_getrfnp
  31. #undef mkl_dc_lapacke_getrfnp_convert
  32. #undef mkl_dc_getrfnp_1x1
  33. #undef mkl_dc_getrfnp_2x2
  34. #undef mkl_dc_getrfnp_3x3
  35. #undef mkl_dc_getrfnp_4x4
  36. #undef mkl_dc_getrfnp_5x5
  37. #undef mkl_dc_getrfnp_mxn
  38. #undef mkl_dc_getrfnp_generic
  39. #undef mkl_dc_getri
  40. #undef mkl_dc_lapacke_getri_convert
  41. #undef mkl_dc_getri_1x1
  42. #undef mkl_dc_getri_2x2
  43. #undef mkl_dc_getri_3x3
  44. #undef mkl_dc_getri_4x4
  45. #undef mkl_dc_getri_5x5
  46. #undef mkl_dc_getri_6x6
  47. #undef mkl_dc_getri_n
  48. #undef mkl_dc_getri_inverse_upper
  49. #undef mkl_dc_getri_solve_lower
  50. #undef mkl_dc_getri_generic
  51. #undef mkl_dc_getrs
  52. #undef mkl_dc_lapacke_getrs_convert
  53. #undef mkl_dc_getrs_trans
  54. #undef mkl_dc_getrs_nxnrhs
  55. #undef mkl_dc_getrs_generic
  56. #undef mkl_dc_geqrf
  57. #undef mkl_dc_lapacke_geqrf_convert
  58. #undef mkl_dc_geqrf_mxn
  59. #undef mkl_dc_geqrf_generic
  60. #undef mkl_dc_geqrf_hh
  61. #undef mkl_dc_geqrf_apply
  62. #undef mkl_dc_potrf
  63. #undef mkl_dc_lapacke_potrf_convert
  64. #undef mkl_dc_potrf_1
  65. #undef mkl_dc_potrf_2
  66. #undef mkl_dc_potrf_3
  67. #undef mkl_dc_potrf_4
  68. #undef mkl_dc_potrf_5
  69. #undef mkl_dc_potrf_n
  70. #undef mkl_dc_potrf_generic
  71. #undef mkl_dc_potrf_uplo
  72. #if defined(MKL_DOUBLE)
  73. #define mkl_dc_getrf mkl_dc_dgetrf
  74. #define mkl_dc_lapacke_getrf_convert mkl_dc_lapacke_dgetrf_convert
  75. #define mkl_dc_getrfnp mkl_dc_dgetrfnp
  76. #define mkl_dc_lapacke_getrfnp_convert mkl_dc_lapacke_dgetrfnp_convert
  77. #define mkl_dc_getri mkl_dc_dgetri
  78. #define mkl_dc_lapacke_getri_convert mkl_dc_lapacke_dgetri_convert
  79. #define mkl_dc_getrs mkl_dc_dgetrs
  80. #define mkl_dc_lapacke_getrs_convert mkl_dc_lapacke_dgetrs_convert
  81. #define mkl_dc_geqrf mkl_dc_dgeqrf
  82. #define mkl_dc_lapacke_geqrf_convert mkl_dc_lapacke_dgeqrf_convert
  83. #define mkl_dc_potrf mkl_dc_dpotrf
  84. #define mkl_dc_lapacke_potrf_convert mkl_dc_lapacke_dpotrf_convert
  85. #elif defined(MKL_SINGLE)
  86. #define mkl_dc_getrf mkl_dc_sgetrf
  87. #define mkl_dc_lapacke_getrf_convert mkl_dc_lapacke_sgetrf_convert
  88. #define mkl_dc_getrfnp mkl_dc_sgetrfnp
  89. #define mkl_dc_lapacke_getrfnp_convert mkl_dc_lapacke_sgetrfnp_convert
  90. #define mkl_dc_getri mkl_dc_sgetri
  91. #define mkl_dc_lapacke_getri_convert mkl_dc_lapacke_sgetri_convert
  92. #define mkl_dc_getrs mkl_dc_sgetrs
  93. #define mkl_dc_lapacke_getrs_convert mkl_dc_lapacke_sgetrs_convert
  94. #define mkl_dc_geqrf mkl_dc_sgeqrf
  95. #define mkl_dc_lapacke_geqrf_convert mkl_dc_lapacke_sgeqrf_convert
  96. #define mkl_dc_potrf mkl_dc_spotrf
  97. #define mkl_dc_lapacke_potrf_convert mkl_dc_lapacke_spotrf_convert
  98. #elif defined(MKL_COMPLEX)
  99. #define mkl_dc_getrf mkl_dc_cgetrf
  100. #define mkl_dc_lapacke_getrf_convert mkl_dc_lapacke_cgetrf_convert
  101. #define mkl_dc_getrfnp mkl_dc_cgetrfnp
  102. #define mkl_dc_lapacke_getrfnp_convert mkl_dc_lapacke_cgetrfnp_convert
  103. #define mkl_dc_getri mkl_dc_cgetri
  104. #define mkl_dc_lapacke_getri_convert mkl_dc_lapacke_cgetri_convert
  105. #define mkl_dc_getrs mkl_dc_cgetrs
  106. #define mkl_dc_lapacke_getrs_convert mkl_dc_lapacke_cgetrs_convert
  107. #define mkl_dc_geqrf mkl_dc_cgeqrf
  108. #define mkl_dc_lapacke_geqrf_convert mkl_dc_lapacke_cgeqrf_convert
  109. #define mkl_dc_potrf mkl_dc_cpotrf
  110. #define mkl_dc_lapacke_potrf_convert mkl_dc_lapacke_cpotrf_convert
  111. #elif defined(MKL_COMPLEX16)
  112. #define mkl_dc_getrf mkl_dc_zgetrf
  113. #define mkl_dc_lapacke_getrf_convert mkl_dc_lapacke_zgetrf_convert
  114. #define mkl_dc_getrfnp mkl_dc_zgetrfnp
  115. #define mkl_dc_lapacke_getrfnp_convert mkl_dc_lapacke_zgetrfnp_convert
  116. #define mkl_dc_getri mkl_dc_zgetri
  117. #define mkl_dc_lapacke_getri_convert mkl_dc_lapacke_zgetri_convert
  118. #define mkl_dc_getrs mkl_dc_zgetrs
  119. #define mkl_dc_lapacke_getrs_convert mkl_dc_lapacke_zgetrs_convert
  120. #define mkl_dc_geqrf mkl_dc_zgeqrf
  121. #define mkl_dc_lapacke_geqrf_convert mkl_dc_lapacke_zgeqrf_convert
  122. #define mkl_dc_potrf mkl_dc_zpotrf
  123. #define mkl_dc_lapacke_potrf_convert mkl_dc_lapacke_zpotrf_convert
  124. #endif
  125. /* ?GETRF */
  126. #define mkl_dc_getrf_mxn(m, n, a, lda, ipiv, info, dc_access) \
  127. do { \
  128. MKL_INT i, j, k, ii, jj; \
  129. *info = 0; \
  130. for (j = 0; j < MKL_DC_MIN(m, n); ++j) { \
  131. mkl_dc_real_type pivot = MKL_DC_ABS1(dc_access(a, lda, j, j)); \
  132. MKL_INT jp = j; \
  133. for (i = j + 1; i < m; ++i) { \
  134. mkl_dc_real_type d = MKL_DC_ABS1(dc_access(a, lda, i, j)); \
  135. if (d > pivot) { \
  136. pivot = d; \
  137. jp = i; \
  138. } \
  139. } \
  140. ipiv[j] = jp + 1; \
  141. if (!MKL_DC_IS_ZERO(dc_access(a, lda, jp, j))) { \
  142. if (jp != j) { \
  143. for (k = 0; k < n; k++) { \
  144. MKL_DC_SWAP(dc_access(a, lda, j, k), dc_access(a, lda, jp, k)); \
  145. } \
  146. } \
  147. for (i = j + 1; i < m; ++i) { \
  148. MKL_DC_DIV(dc_access(a, lda, i, j), dc_access(a, lda, i, j), dc_access(a, lda, j, j)); \
  149. } \
  150. } else if (*info == 0) { \
  151. *info = j + 1; \
  152. } \
  153. for (jj = j + 1; jj < n; ++jj) { \
  154. mkl_dc_type a_j_jj = dc_access(a, lda, j, jj); \
  155. MKL_DC_PRAGMA_VECTOR \
  156. for (ii = j + 1; ii < m; ++ii) { \
  157. MKL_DC_SUB_MUL(dc_access(a, lda, ii, jj), dc_access(a, lda, ii, jj), \
  158. dc_access(a, lda, ii, j), a_j_jj); \
  159. } \
  160. } \
  161. } \
  162. } while (0) \
  163. #define mkl_dc_getrf_1x1(m, n, a, lda, ipiv, info, dc_access) \
  164. do { \
  165. ipiv[0] = 1; \
  166. *info = !MKL_DC_IS_ZERO(dc_access(a, lda, 0, 0)) ? 0 : 1; \
  167. } while (0) \
  168. #define mkl_dc_getrf_2x2(m, n, a, lda, ipiv, info, dc_access) \
  169. do { \
  170. mkl_dc_type a00 = dc_access(a, lda, 0, 0); \
  171. mkl_dc_type a10 = dc_access(a, lda, 1, 0); \
  172. mkl_dc_type a01 = dc_access(a, lda, 0, 1); \
  173. mkl_dc_type a11 = dc_access(a, lda, 1, 1); \
  174. mkl_dc_real_type pivot = MKL_DC_ABS1(a00); \
  175. mkl_dc_real_type d = MKL_DC_ABS1(a10); \
  176. MKL_INT jp = 0; \
  177. *info = 0; \
  178. if (d > pivot) { \
  179. pivot = d; \
  180. jp = 1; \
  181. } \
  182. ipiv[0] = jp + 1; \
  183. if (jp == 0) { \
  184. if (!MKL_DC_IS_ZERO(a00)) { \
  185. MKL_DC_DIV(a10, a10, a00); \
  186. MKL_DC_SUB_MUL(a11, a11, a10, a01); \
  187. } else { \
  188. *info = 1; \
  189. } \
  190. } else { \
  191. MKL_DC_SWAP(a00, a10); \
  192. MKL_DC_SWAP(a01, a11); \
  193. MKL_DC_DIV(a10, a10, a00); \
  194. MKL_DC_SUB_MUL(a11, a11, a10, a01); \
  195. } \
  196. ipiv[1] = 2; \
  197. if (MKL_DC_IS_ZERO(a11) && *info == 0) { \
  198. *info = 2; \
  199. } \
  200. dc_access(a, lda, 0, 0) = a00; \
  201. dc_access(a, lda, 1, 0) = a10; \
  202. dc_access(a, lda, 0, 1) = a01; \
  203. dc_access(a, lda, 1, 1) = a11; \
  204. } while (0) \
  205. #define mkl_dc_getrf_3x3(m, n, a, lda, ipiv, info, dc_access) \
  206. do { \
  207. mkl_dc_type a00 = dc_access(a, lda, 0, 0); \
  208. mkl_dc_type a10 = dc_access(a, lda, 1, 0); \
  209. mkl_dc_type a20 = dc_access(a, lda, 2, 0); \
  210. mkl_dc_type a01 = dc_access(a, lda, 0, 1); \
  211. mkl_dc_type a11 = dc_access(a, lda, 1, 1); \
  212. mkl_dc_type a21 = dc_access(a, lda, 2, 1); \
  213. mkl_dc_type a02 = dc_access(a, lda, 0, 2); \
  214. mkl_dc_type a12 = dc_access(a, lda, 1, 2); \
  215. mkl_dc_type a22 = dc_access(a, lda, 2, 2); \
  216. mkl_dc_real_type pivot = MKL_DC_ABS1(a00); \
  217. mkl_dc_real_type d = MKL_DC_ABS1(a10); \
  218. MKL_INT jp = 0; \
  219. *info = 0; \
  220. if (d > pivot) { \
  221. pivot = d; \
  222. jp = 1; \
  223. } \
  224. d = MKL_DC_ABS1(a20); \
  225. if (d > pivot) { \
  226. pivot = d; \
  227. jp = 2; \
  228. } \
  229. ipiv[0] = jp + 1; \
  230. switch (jp) { \
  231. case 0: \
  232. if (!MKL_DC_IS_ZERO(a00)) { \
  233. MKL_DC_DIV(a10, a10, a00); \
  234. MKL_DC_DIV(a20, a20, a00); \
  235. } else { \
  236. *info = 1; \
  237. } \
  238. break; \
  239. case 1: \
  240. if (!MKL_DC_IS_ZERO(a10)) { \
  241. MKL_DC_SWAP(a00, a10); \
  242. MKL_DC_SWAP(a01, a11); \
  243. MKL_DC_SWAP(a02, a12); \
  244. MKL_DC_DIV(a10, a10, a00); \
  245. MKL_DC_DIV(a20, a20, a00); \
  246. } else { \
  247. *info = 1; \
  248. } \
  249. break; \
  250. case 2: \
  251. if (!MKL_DC_IS_ZERO(a20)) { \
  252. MKL_DC_SWAP(a00, a20); \
  253. MKL_DC_SWAP(a01, a21); \
  254. MKL_DC_SWAP(a02, a22); \
  255. MKL_DC_DIV(a10, a10, a00); \
  256. MKL_DC_DIV(a20, a20, a00); \
  257. } else { \
  258. *info = 1; \
  259. } \
  260. break; \
  261. } \
  262. MKL_DC_SUB_MUL(a11, a11, a10, a01); \
  263. MKL_DC_SUB_MUL(a21, a21, a20, a01); \
  264. MKL_DC_SUB_MUL(a12, a12, a10, a02); \
  265. MKL_DC_SUB_MUL(a22, a22, a20, a02); \
  266. pivot = MKL_DC_ABS1(a11); \
  267. jp = 1; \
  268. d = MKL_DC_ABS1(a21); \
  269. if (d > pivot) { \
  270. pivot = d; \
  271. jp = 2; \
  272. } \
  273. ipiv[1] = jp + 1; \
  274. switch (jp) { \
  275. case 1: \
  276. if (!MKL_DC_IS_ZERO(a11)) { \
  277. MKL_DC_DIV(a21, a21, a11); \
  278. } else if (*info == 0) { \
  279. *info = 2; \
  280. } \
  281. break; \
  282. case 2: \
  283. if (!MKL_DC_IS_ZERO(a21)) { \
  284. MKL_DC_SWAP(a10, a20); \
  285. MKL_DC_SWAP(a11, a21); \
  286. MKL_DC_SWAP(a12, a22); \
  287. MKL_DC_DIV(a21, a21, a11); \
  288. } else if (*info == 0) { \
  289. *info = 2; \
  290. } \
  291. break; \
  292. } \
  293. MKL_DC_SUB_MUL(a22, a22, a21, a12); \
  294. ipiv[2] = 3; \
  295. if (MKL_DC_IS_ZERO(a22) && *info == 0) { \
  296. *info = 3; \
  297. } \
  298. dc_access(a, lda, 0, 0) = a00; \
  299. dc_access(a, lda, 1, 0) = a10; \
  300. dc_access(a, lda, 2, 0) = a20; \
  301. dc_access(a, lda, 0, 1) = a01; \
  302. dc_access(a, lda, 1, 1) = a11; \
  303. dc_access(a, lda, 2, 1) = a21; \
  304. dc_access(a, lda, 0, 2) = a02; \
  305. dc_access(a, lda, 1, 2) = a12; \
  306. dc_access(a, lda, 2, 2) = a22; \
  307. } while (0) \
  308. #define mkl_dc_getrf_4x4(m, n, a, lda, ipiv, info, dc_access) \
  309. do { \
  310. mkl_dc_type a00 = dc_access(a, lda, 0, 0); \
  311. mkl_dc_type a10 = dc_access(a, lda, 1, 0); \
  312. mkl_dc_type a20 = dc_access(a, lda, 2, 0); \
  313. mkl_dc_type a30 = dc_access(a, lda, 3, 0); \
  314. mkl_dc_type a01 = dc_access(a, lda, 0, 1); \
  315. mkl_dc_type a11 = dc_access(a, lda, 1, 1); \
  316. mkl_dc_type a21 = dc_access(a, lda, 2, 1); \
  317. mkl_dc_type a31 = dc_access(a, lda, 3, 1); \
  318. mkl_dc_type a02 = dc_access(a, lda, 0, 2); \
  319. mkl_dc_type a12 = dc_access(a, lda, 1, 2); \
  320. mkl_dc_type a22 = dc_access(a, lda, 2, 2); \
  321. mkl_dc_type a32 = dc_access(a, lda, 3, 2); \
  322. mkl_dc_type a03 = dc_access(a, lda, 0, 3); \
  323. mkl_dc_type a13 = dc_access(a, lda, 1, 3); \
  324. mkl_dc_type a23 = dc_access(a, lda, 2, 3); \
  325. mkl_dc_type a33 = dc_access(a, lda, 3, 3); \
  326. mkl_dc_real_type pivot = MKL_DC_ABS1(a00); \
  327. mkl_dc_real_type d = MKL_DC_ABS1(a10); \
  328. MKL_INT jp = 0; \
  329. *info = 0; \
  330. if (d > pivot) { \
  331. pivot = d; \
  332. jp = 1; \
  333. } \
  334. d = MKL_DC_ABS1(a20); \
  335. if (d > pivot) { \
  336. pivot = d; \
  337. jp = 2; \
  338. } \
  339. d = MKL_DC_ABS1(a30); \
  340. if (d > pivot) { \
  341. pivot = d; \
  342. jp = 3; \
  343. } \
  344. ipiv[0] = jp + 1; \
  345. switch (jp) { \
  346. case 0: \
  347. if (!MKL_DC_IS_ZERO(a00)) { \
  348. MKL_DC_DIV(a10, a10, a00); \
  349. MKL_DC_DIV(a20, a20, a00); \
  350. MKL_DC_DIV(a30, a30, a00); \
  351. } else { \
  352. *info = 1; \
  353. } \
  354. break; \
  355. case 1: \
  356. if (!MKL_DC_IS_ZERO(a10)) { \
  357. MKL_DC_SWAP(a00, a10); \
  358. MKL_DC_SWAP(a01, a11); \
  359. MKL_DC_SWAP(a02, a12); \
  360. MKL_DC_SWAP(a03, a13); \
  361. MKL_DC_DIV(a10, a10, a00); \
  362. MKL_DC_DIV(a20, a20, a00); \
  363. MKL_DC_DIV(a30, a30, a00); \
  364. } else { \
  365. *info = 1; \
  366. } \
  367. break; \
  368. case 2: \
  369. if (!MKL_DC_IS_ZERO(a20)) { \
  370. MKL_DC_SWAP(a00, a20); \
  371. MKL_DC_SWAP(a01, a21); \
  372. MKL_DC_SWAP(a02, a22); \
  373. MKL_DC_SWAP(a03, a23); \
  374. MKL_DC_DIV(a10, a10, a00); \
  375. MKL_DC_DIV(a20, a20, a00); \
  376. MKL_DC_DIV(a30, a30, a00); \
  377. } else { \
  378. *info = 1; \
  379. } \
  380. break; \
  381. case 3: \
  382. if (!MKL_DC_IS_ZERO(a30)) { \
  383. MKL_DC_SWAP(a00, a30); \
  384. MKL_DC_SWAP(a01, a31); \
  385. MKL_DC_SWAP(a02, a32); \
  386. MKL_DC_SWAP(a03, a33); \
  387. MKL_DC_DIV(a10, a10, a00); \
  388. MKL_DC_DIV(a20, a20, a00); \
  389. MKL_DC_DIV(a30, a30, a00); \
  390. } else { \
  391. *info = 1; \
  392. } \
  393. break; \
  394. } \
  395. MKL_DC_SUB_MUL(a11, a11, a10, a01); \
  396. MKL_DC_SUB_MUL(a21, a21, a20, a01); \
  397. MKL_DC_SUB_MUL(a31, a31, a30, a01); \
  398. MKL_DC_SUB_MUL(a12, a12, a10, a02); \
  399. MKL_DC_SUB_MUL(a22, a22, a20, a02); \
  400. MKL_DC_SUB_MUL(a32, a32, a30, a02); \
  401. MKL_DC_SUB_MUL(a13, a13, a10, a03); \
  402. MKL_DC_SUB_MUL(a23, a23, a20, a03); \
  403. MKL_DC_SUB_MUL(a33, a33, a30, a03); \
  404. pivot = MKL_DC_ABS1(a11); \
  405. jp = 1; \
  406. d = MKL_DC_ABS1(a21); \
  407. if (d > pivot) { \
  408. pivot = d; \
  409. jp = 2; \
  410. } \
  411. d = MKL_DC_ABS1(a31); \
  412. if (d > pivot) { \
  413. pivot = d; \
  414. jp = 3; \
  415. } \
  416. ipiv[1] = jp + 1; \
  417. switch (jp) { \
  418. case 1: \
  419. if (!MKL_DC_IS_ZERO(a11)) { \
  420. MKL_DC_DIV(a21, a21, a11); \
  421. MKL_DC_DIV(a31, a31, a11); \
  422. } else if (*info == 0) { \
  423. *info = 2; \
  424. } \
  425. break; \
  426. case 2: \
  427. if (!MKL_DC_IS_ZERO(a21)) { \
  428. MKL_DC_SWAP(a10, a20); \
  429. MKL_DC_SWAP(a11, a21); \
  430. MKL_DC_SWAP(a12, a22); \
  431. MKL_DC_SWAP(a13, a23); \
  432. MKL_DC_DIV(a21, a21, a11); \
  433. MKL_DC_DIV(a31, a31, a11); \
  434. } else if (*info == 0) { \
  435. *info = 2; \
  436. } \
  437. break; \
  438. case 3: \
  439. if (!MKL_DC_IS_ZERO(a31)) { \
  440. MKL_DC_SWAP(a10, a30); \
  441. MKL_DC_SWAP(a11, a31); \
  442. MKL_DC_SWAP(a12, a32); \
  443. MKL_DC_SWAP(a13, a33); \
  444. MKL_DC_DIV(a21, a21, a11); \
  445. MKL_DC_DIV(a31, a31, a11); \
  446. } else if (*info == 0) { \
  447. *info = 2; \
  448. } \
  449. break; \
  450. } \
  451. MKL_DC_SUB_MUL(a22, a22, a21, a12); \
  452. MKL_DC_SUB_MUL(a32, a32, a31, a12); \
  453. MKL_DC_SUB_MUL(a23, a23, a21, a13); \
  454. MKL_DC_SUB_MUL(a33, a33, a31, a13); \
  455. pivot = MKL_DC_ABS1(a22); \
  456. jp = 2; \
  457. d = MKL_DC_ABS1(a32); \
  458. if (d > pivot) { \
  459. pivot = d; \
  460. jp = 3; \
  461. } \
  462. ipiv[2] = jp + 1; \
  463. switch (jp) { \
  464. case 2: \
  465. if (!MKL_DC_IS_ZERO(a22)) { \
  466. MKL_DC_DIV(a32, a32, a22); \
  467. } else if (*info == 0) { \
  468. *info = 3; \
  469. } \
  470. break; \
  471. case 3: \
  472. if (!MKL_DC_IS_ZERO(a32)) { \
  473. MKL_DC_SWAP(a20, a30); \
  474. MKL_DC_SWAP(a21, a31); \
  475. MKL_DC_SWAP(a22, a32); \
  476. MKL_DC_SWAP(a23, a33); \
  477. MKL_DC_DIV(a32, a32, a22); \
  478. } else if (*info == 0) { \
  479. *info = 3; \
  480. } \
  481. break; \
  482. } \
  483. MKL_DC_SUB_MUL(a33, a33, a32, a23); \
  484. ipiv[3] = 4; \
  485. if (MKL_DC_IS_ZERO(a33) && *info == 0) { \
  486. *info = 4; \
  487. } \
  488. dc_access(a, lda, 0, 0) = a00; \
  489. dc_access(a, lda, 1, 0) = a10; \
  490. dc_access(a, lda, 2, 0) = a20; \
  491. dc_access(a, lda, 3, 0) = a30; \
  492. dc_access(a, lda, 0, 1) = a01; \
  493. dc_access(a, lda, 1, 1) = a11; \
  494. dc_access(a, lda, 2, 1) = a21; \
  495. dc_access(a, lda, 3, 1) = a31; \
  496. dc_access(a, lda, 0, 2) = a02; \
  497. dc_access(a, lda, 1, 2) = a12; \
  498. dc_access(a, lda, 2, 2) = a22; \
  499. dc_access(a, lda, 3, 2) = a32; \
  500. dc_access(a, lda, 0, 3) = a03; \
  501. dc_access(a, lda, 1, 3) = a13; \
  502. dc_access(a, lda, 2, 3) = a23; \
  503. dc_access(a, lda, 3, 3) = a33; \
  504. } while (0)
  505. #define mkl_dc_getrf_5x5(m, n, a, lda, ipiv, info, dc_access) \
  506. do { \
  507. mkl_dc_type a00 = dc_access(a, lda, 0, 0); \
  508. mkl_dc_type a10 = dc_access(a, lda, 1, 0); \
  509. mkl_dc_type a20 = dc_access(a, lda, 2, 0); \
  510. mkl_dc_type a30 = dc_access(a, lda, 3, 0); \
  511. mkl_dc_type a40 = dc_access(a, lda, 4, 0); \
  512. mkl_dc_type a01 = dc_access(a, lda, 0, 1); \
  513. mkl_dc_type a11 = dc_access(a, lda, 1, 1); \
  514. mkl_dc_type a21 = dc_access(a, lda, 2, 1); \
  515. mkl_dc_type a31 = dc_access(a, lda, 3, 1); \
  516. mkl_dc_type a41 = dc_access(a, lda, 4, 1); \
  517. mkl_dc_type a02 = dc_access(a, lda, 0, 2); \
  518. mkl_dc_type a12 = dc_access(a, lda, 1, 2); \
  519. mkl_dc_type a22 = dc_access(a, lda, 2, 2); \
  520. mkl_dc_type a32 = dc_access(a, lda, 3, 2); \
  521. mkl_dc_type a42 = dc_access(a, lda, 4, 2); \
  522. mkl_dc_type a03 = dc_access(a, lda, 0, 3); \
  523. mkl_dc_type a13 = dc_access(a, lda, 1, 3); \
  524. mkl_dc_type a23 = dc_access(a, lda, 2, 3); \
  525. mkl_dc_type a33 = dc_access(a, lda, 3, 3); \
  526. mkl_dc_type a43 = dc_access(a, lda, 4, 3); \
  527. mkl_dc_type a04 = dc_access(a, lda, 0, 4); \
  528. mkl_dc_type a14 = dc_access(a, lda, 1, 4); \
  529. mkl_dc_type a24 = dc_access(a, lda, 2, 4); \
  530. mkl_dc_type a34 = dc_access(a, lda, 3, 4); \
  531. mkl_dc_type a44 = dc_access(a, lda, 4, 4); \
  532. mkl_dc_real_type pivot = MKL_DC_ABS1(a00); \
  533. mkl_dc_real_type d = MKL_DC_ABS1(a10); \
  534. MKL_INT jp = 0; \
  535. *info = 0; \
  536. if (d > pivot) { \
  537. pivot = d; \
  538. jp = 1; \
  539. } \
  540. d = MKL_DC_ABS1(a20); \
  541. if (d > pivot) { \
  542. pivot = d; \
  543. jp = 2; \
  544. } \
  545. d = MKL_DC_ABS1(a30); \
  546. if (d > pivot) { \
  547. pivot = d; \
  548. jp = 3; \
  549. } \
  550. d = MKL_DC_ABS1(a40); \
  551. if (d > pivot) { \
  552. pivot = d; \
  553. jp = 4; \
  554. } \
  555. ipiv[0] = jp + 1; \
  556. switch (jp) { \
  557. case 0: \
  558. if (!MKL_DC_IS_ZERO(a00)) { \
  559. MKL_DC_DIV(a10, a10, a00); \
  560. MKL_DC_DIV(a20, a20, a00); \
  561. MKL_DC_DIV(a30, a30, a00); \
  562. MKL_DC_DIV(a40, a40, a00); \
  563. } else { \
  564. *info = 1; \
  565. } \
  566. break; \
  567. case 1: \
  568. if (!MKL_DC_IS_ZERO(a10)) { \
  569. MKL_DC_SWAP(a00, a10); \
  570. MKL_DC_SWAP(a01, a11); \
  571. MKL_DC_SWAP(a02, a12); \
  572. MKL_DC_SWAP(a03, a13); \
  573. MKL_DC_SWAP(a04, a14); \
  574. MKL_DC_DIV(a10, a10, a00); \
  575. MKL_DC_DIV(a20, a20, a00); \
  576. MKL_DC_DIV(a30, a30, a00); \
  577. MKL_DC_DIV(a40, a40, a00); \
  578. } else { \
  579. *info = 1; \
  580. } \
  581. break; \
  582. case 2: \
  583. if (!MKL_DC_IS_ZERO(a20)) { \
  584. MKL_DC_SWAP(a00, a20); \
  585. MKL_DC_SWAP(a01, a21); \
  586. MKL_DC_SWAP(a02, a22); \
  587. MKL_DC_SWAP(a03, a23); \
  588. MKL_DC_SWAP(a04, a24); \
  589. MKL_DC_DIV(a10, a10, a00); \
  590. MKL_DC_DIV(a20, a20, a00); \
  591. MKL_DC_DIV(a30, a30, a00); \
  592. MKL_DC_DIV(a40, a40, a00); \
  593. } else { \
  594. *info = 1; \
  595. } \
  596. break; \
  597. case 3: \
  598. if (!MKL_DC_IS_ZERO(a30)) { \
  599. MKL_DC_SWAP(a00, a30); \
  600. MKL_DC_SWAP(a01, a31); \
  601. MKL_DC_SWAP(a02, a32); \
  602. MKL_DC_SWAP(a03, a33); \
  603. MKL_DC_SWAP(a04, a34); \
  604. MKL_DC_DIV(a10, a10, a00); \
  605. MKL_DC_DIV(a20, a20, a00); \
  606. MKL_DC_DIV(a30, a30, a00); \
  607. MKL_DC_DIV(a40, a40, a00); \
  608. } else { \
  609. *info = 1; \
  610. } \
  611. break; \
  612. case 4: \
  613. if (!MKL_DC_IS_ZERO(a40)) { \
  614. MKL_DC_SWAP(a00, a40); \
  615. MKL_DC_SWAP(a01, a41); \
  616. MKL_DC_SWAP(a02, a42); \
  617. MKL_DC_SWAP(a03, a43); \
  618. MKL_DC_SWAP(a04, a44); \
  619. MKL_DC_DIV(a10, a10, a00); \
  620. MKL_DC_DIV(a20, a20, a00); \
  621. MKL_DC_DIV(a30, a30, a00); \
  622. MKL_DC_DIV(a40, a40, a00); \
  623. } else { \
  624. *info = 1; \
  625. } \
  626. break; \
  627. } \
  628. MKL_DC_SUB_MUL(a11, a11, a10, a01); \
  629. MKL_DC_SUB_MUL(a21, a21, a20, a01); \
  630. MKL_DC_SUB_MUL(a31, a31, a30, a01); \
  631. MKL_DC_SUB_MUL(a41, a41, a40, a01); \
  632. MKL_DC_SUB_MUL(a12, a12, a10, a02); \
  633. MKL_DC_SUB_MUL(a22, a22, a20, a02); \
  634. MKL_DC_SUB_MUL(a32, a32, a30, a02); \
  635. MKL_DC_SUB_MUL(a42, a42, a40, a02); \
  636. MKL_DC_SUB_MUL(a13, a13, a10, a03); \
  637. MKL_DC_SUB_MUL(a23, a23, a20, a03); \
  638. MKL_DC_SUB_MUL(a33, a33, a30, a03); \
  639. MKL_DC_SUB_MUL(a43, a43, a40, a03); \
  640. MKL_DC_SUB_MUL(a14, a14, a10, a04); \
  641. MKL_DC_SUB_MUL(a24, a24, a20, a04); \
  642. MKL_DC_SUB_MUL(a34, a34, a30, a04); \
  643. MKL_DC_SUB_MUL(a44, a44, a40, a04); \
  644. pivot = MKL_DC_ABS1(a11); \
  645. jp = 1; \
  646. d = MKL_DC_ABS1(a21); \
  647. if (d > pivot) { \
  648. pivot = d; \
  649. jp = 2; \
  650. } \
  651. d = MKL_DC_ABS1(a31); \
  652. if (d > pivot) { \
  653. pivot = d; \
  654. jp = 3; \
  655. } \
  656. d = MKL_DC_ABS1(a41); \
  657. if (d > pivot) { \
  658. pivot = d; \
  659. jp = 4; \
  660. } \
  661. ipiv[1] = jp + 1; \
  662. switch (jp) { \
  663. case 1: \
  664. if (!MKL_DC_IS_ZERO(a11)) { \
  665. MKL_DC_DIV(a21, a21, a11); \
  666. MKL_DC_DIV(a31, a31, a11); \
  667. MKL_DC_DIV(a41, a41, a11); \
  668. } else if (*info == 0) { \
  669. *info = 2; \
  670. } \
  671. break; \
  672. case 2: \
  673. if (!MKL_DC_IS_ZERO(a21)) { \
  674. MKL_DC_SWAP(a10, a20); \
  675. MKL_DC_SWAP(a11, a21); \
  676. MKL_DC_SWAP(a12, a22); \
  677. MKL_DC_SWAP(a13, a23); \
  678. MKL_DC_SWAP(a14, a24); \
  679. MKL_DC_DIV(a21, a21, a11); \
  680. MKL_DC_DIV(a31, a31, a11); \
  681. MKL_DC_DIV(a41, a41, a11); \
  682. } else if (*info == 0) { \
  683. *info = 2; \
  684. } \
  685. break; \
  686. case 3: \
  687. if (!MKL_DC_IS_ZERO(a31)) { \
  688. MKL_DC_SWAP(a10, a30); \
  689. MKL_DC_SWAP(a11, a31); \
  690. MKL_DC_SWAP(a12, a32); \
  691. MKL_DC_SWAP(a13, a33); \
  692. MKL_DC_SWAP(a14, a34); \
  693. MKL_DC_DIV(a21, a21, a11); \
  694. MKL_DC_DIV(a31, a31, a11); \
  695. MKL_DC_DIV(a41, a41, a11); \
  696. } else if (*info == 0) { \
  697. *info = 2; \
  698. } \
  699. break; \
  700. case 4: \
  701. if (!MKL_DC_IS_ZERO(a41)) { \
  702. MKL_DC_SWAP(a10, a40); \
  703. MKL_DC_SWAP(a11, a41); \
  704. MKL_DC_SWAP(a12, a42); \
  705. MKL_DC_SWAP(a13, a43); \
  706. MKL_DC_SWAP(a14, a44); \
  707. MKL_DC_DIV(a21, a21, a11); \
  708. MKL_DC_DIV(a31, a31, a11); \
  709. MKL_DC_DIV(a41, a41, a11); \
  710. } else if (*info == 0) { \
  711. *info = 2; \
  712. } \
  713. break; \
  714. } \
  715. MKL_DC_SUB_MUL(a22, a22, a21, a12); \
  716. MKL_DC_SUB_MUL(a32, a32, a31, a12); \
  717. MKL_DC_SUB_MUL(a42, a42, a41, a12); \
  718. MKL_DC_SUB_MUL(a23, a23, a21, a13); \
  719. MKL_DC_SUB_MUL(a33, a33, a31, a13); \
  720. MKL_DC_SUB_MUL(a43, a43, a41, a13); \
  721. MKL_DC_SUB_MUL(a24, a24, a21, a14); \
  722. MKL_DC_SUB_MUL(a34, a34, a31, a14); \
  723. MKL_DC_SUB_MUL(a44, a44, a41, a14); \
  724. pivot = MKL_DC_ABS1(a22); \
  725. jp = 2; \
  726. d = MKL_DC_ABS1(a32); \
  727. if (d > pivot) { \
  728. pivot = d; \
  729. jp = 3; \
  730. } \
  731. d = MKL_DC_ABS1(a42); \
  732. if (d > pivot) { \
  733. pivot = d; \
  734. jp = 4; \
  735. } \
  736. ipiv[2] = jp + 1; \
  737. switch (jp) { \
  738. case 2: \
  739. if (!MKL_DC_IS_ZERO(a22)) { \
  740. MKL_DC_DIV(a32, a32, a22); \
  741. MKL_DC_DIV(a42, a42, a22); \
  742. } else if (*info == 0) { \
  743. *info = 3; \
  744. } \
  745. break; \
  746. case 3: \
  747. if (!MKL_DC_IS_ZERO(a32)) { \
  748. MKL_DC_SWAP(a20, a30); \
  749. MKL_DC_SWAP(a21, a31); \
  750. MKL_DC_SWAP(a22, a32); \
  751. MKL_DC_SWAP(a23, a33); \
  752. MKL_DC_SWAP(a24, a34); \
  753. MKL_DC_DIV(a32, a32, a22); \
  754. MKL_DC_DIV(a42, a42, a22); \
  755. } else if (*info == 0) { \
  756. *info = 3; \
  757. } \
  758. break; \
  759. case 4: \
  760. if (!MKL_DC_IS_ZERO(a42)) { \
  761. MKL_DC_SWAP(a20, a40); \
  762. MKL_DC_SWAP(a21, a41); \
  763. MKL_DC_SWAP(a22, a42); \
  764. MKL_DC_SWAP(a23, a43); \
  765. MKL_DC_SWAP(a24, a44); \
  766. MKL_DC_DIV(a32, a32, a22); \
  767. MKL_DC_DIV(a42, a42, a22); \
  768. } else if (*info == 0) { \
  769. *info = 3; \
  770. } \
  771. break; \
  772. } \
  773. MKL_DC_SUB_MUL(a33, a33, a32, a23); \
  774. MKL_DC_SUB_MUL(a43, a43, a42, a23); \
  775. MKL_DC_SUB_MUL(a34, a34, a32, a24); \
  776. MKL_DC_SUB_MUL(a44, a44, a42, a24); \
  777. pivot = MKL_DC_ABS1(a33); \
  778. jp = 3; \
  779. d = MKL_DC_ABS1(a43); \
  780. if (d > pivot) { \
  781. pivot = d; \
  782. jp = 4; \
  783. } \
  784. ipiv[3] = jp + 1; \
  785. switch (jp) { \
  786. case 3: \
  787. if (!MKL_DC_IS_ZERO(a33)) { \
  788. MKL_DC_DIV(a43, a43, a33); \
  789. } else if (*info == 0) { \
  790. *info = 4; \
  791. } \
  792. break; \
  793. case 4: \
  794. if (!MKL_DC_IS_ZERO(a43)) { \
  795. MKL_DC_SWAP(a30, a40); \
  796. MKL_DC_SWAP(a31, a41); \
  797. MKL_DC_SWAP(a32, a42); \
  798. MKL_DC_SWAP(a33, a43); \
  799. MKL_DC_SWAP(a34, a44); \
  800. MKL_DC_DIV(a43, a43, a33); \
  801. } else if (*info == 0) { \
  802. *info = 4; \
  803. } \
  804. break; \
  805. } \
  806. MKL_DC_SUB_MUL(a44, a44, a43, a34); \
  807. ipiv[4] = 5; \
  808. if (MKL_DC_IS_ZERO(a44) && *info == 0) { \
  809. *info = 5; \
  810. } \
  811. dc_access(a, lda, 0, 0) = a00; \
  812. dc_access(a, lda, 1, 0) = a10; \
  813. dc_access(a, lda, 2, 0) = a20; \
  814. dc_access(a, lda, 3, 0) = a30; \
  815. dc_access(a, lda, 4, 0) = a40; \
  816. dc_access(a, lda, 0, 1) = a01; \
  817. dc_access(a, lda, 1, 1) = a11; \
  818. dc_access(a, lda, 2, 1) = a21; \
  819. dc_access(a, lda, 3, 1) = a31; \
  820. dc_access(a, lda, 4, 1) = a41; \
  821. dc_access(a, lda, 0, 2) = a02; \
  822. dc_access(a, lda, 1, 2) = a12; \
  823. dc_access(a, lda, 2, 2) = a22; \
  824. dc_access(a, lda, 3, 2) = a32; \
  825. dc_access(a, lda, 4, 2) = a42; \
  826. dc_access(a, lda, 0, 3) = a03; \
  827. dc_access(a, lda, 1, 3) = a13; \
  828. dc_access(a, lda, 2, 3) = a23; \
  829. dc_access(a, lda, 3, 3) = a33; \
  830. dc_access(a, lda, 4, 3) = a43; \
  831. dc_access(a, lda, 0, 4) = a04; \
  832. dc_access(a, lda, 1, 4) = a14; \
  833. dc_access(a, lda, 2, 4) = a24; \
  834. dc_access(a, lda, 3, 4) = a34; \
  835. dc_access(a, lda, 4, 4) = a44; \
  836. } while (0)
  837. #define mkl_dc_getrf_generic(m, n, a, lda, ipiv, info, dc_access) \
  838. do { \
  839. if (m == n && m <= 5) { \
  840. switch (m) { \
  841. case 0: \
  842. break; \
  843. case 1: \
  844. mkl_dc_getrf_1x1(m, n, a, lda, ipiv, info, dc_access); \
  845. break; \
  846. case 2: \
  847. mkl_dc_getrf_2x2(m, n, a, lda, ipiv, info, dc_access); \
  848. break; \
  849. case 3: \
  850. mkl_dc_getrf_3x3(m, n, a, lda, ipiv, info, dc_access); \
  851. break; \
  852. case 4: \
  853. mkl_dc_getrf_4x4(m, n, a, lda, ipiv, info, dc_access); \
  854. break; \
  855. case 5: \
  856. mkl_dc_getrf_5x5(m, n, a, lda, ipiv, info, dc_access); \
  857. break; \
  858. } \
  859. } else { \
  860. mkl_dc_getrf_mxn(m, n, a, lda, ipiv, info, dc_access); \
  861. } \
  862. } while (0)
  863. static __inline void mkl_dc_getrf(MKL_INT m, MKL_INT n, mkl_dc_type* a, MKL_INT lda, MKL_INT* ipiv, MKL_INT* info)
  864. {
  865. mkl_dc_getrf_generic(m, n, a, lda, ipiv, info, MKL_DC_MN);
  866. }
  867. #ifndef MKL_DIRECT_CALL_LAPACKE_DISABLE
  868. static __inline lapack_int mkl_dc_lapacke_getrf_convert(int matrix_layout, lapack_int m, lapack_int n, mkl_dc_type* a, lapack_int lda, lapack_int* ipiv)
  869. {
  870. lapack_int info = 0;
  871. if (MKL_DC_GETRF_CHECKSIZE(m, n)) {
  872. if (matrix_layout == LAPACK_ROW_MAJOR) {
  873. mkl_dc_getrf_generic(m, n, a, lda, ipiv, &info, MKL_DC_MT);
  874. } else {
  875. mkl_dc_getrf_generic(m, n, a, lda, ipiv, &info, MKL_DC_MN);
  876. }
  877. return info;
  878. }
  879. return MKL_DC_CONCAT3(LAPACKE_, MKL_DC_PREC_LETTER, getrf)(matrix_layout, m, n, a, lda, ipiv);
  880. }
  881. #endif
  882. /* ?GETRFNP */
  883. #define mkl_dc_getrfnp_mxn(m, n, a, lda, info, dc_access) \
  884. do { \
  885. MKL_INT i; \
  886. MKL_INT j; \
  887. MKL_INT ii; \
  888. MKL_INT jj; \
  889. MKL_INT is_zero; \
  890. for (j = 0; j < MKL_DC_MIN(m, n); j++) { \
  891. mkl_dc_type ajj = dc_access(a, lda, j, j); \
  892. MKL_DC_IS_ZERO2(is_zero, ajj); \
  893. if (!is_zero) { \
  894. for (i = j + 1; i < m; i++) { \
  895. MKL_DC_DIV(dc_access(a, lda, i, j), dc_access(a, lda, i, j), dc_access(a, lda, j, j)); \
  896. } \
  897. } else if (*info == 0) { \
  898. *info = j + 1; \
  899. } \
  900. for (jj = j + 1; jj < n; jj++) { \
  901. mkl_dc_type a_j_jj = dc_access(a, lda, j, jj); \
  902. for (ii = j + 1; ii < m; ii++) { \
  903. MKL_DC_SUB_MUL(dc_access(a, lda, ii, jj), dc_access(a, lda, ii, jj), dc_access(a, lda, ii, j), a_j_jj); \
  904. } \
  905. } \
  906. } \
  907. } while (0) \
  908. #define mkl_dc_getrfnp_1x1(m, n, a, lda, info, dc_access) \
  909. do { \
  910. MKL_INT is_zero; \
  911. mkl_dc_type a00 = dc_access(a, lda, 0, 0); \
  912. MKL_DC_IS_ZERO2(is_zero, a00); \
  913. *info = !!is_zero; \
  914. } while (0) \
  915. #define mkl_dc_getrfnp_2x2(m, n, a, lda, info, dc_access) \
  916. do { \
  917. MKL_INT is_zero; \
  918. mkl_dc_type a00 = dc_access(a, lda, 0, 0); \
  919. mkl_dc_type a10 = dc_access(a, lda, 1, 0); \
  920. mkl_dc_type a01 = dc_access(a, lda, 0, 1); \
  921. mkl_dc_type a11 = dc_access(a, lda, 1, 1); \
  922. MKL_DC_IS_ZERO2(is_zero, a00); \
  923. if (!is_zero) { \
  924. MKL_DC_DIV(a10, a10, a00); \
  925. MKL_DC_SUB_MUL(a11, a11, a10, a01); \
  926. } else { \
  927. *info = 1; \
  928. } \
  929. MKL_DC_IS_ZERO2(is_zero, a11); \
  930. if (is_zero && *info == 0) { \
  931. *info = 2; \
  932. } \
  933. dc_access(a, lda, 1, 0) = a10; \
  934. dc_access(a, lda, 1, 1) = a11; \
  935. } while (0) \
  936. #define mkl_dc_getrfnp_3x3(m, n, a, lda, info, dc_access) \
  937. do { \
  938. MKL_INT is_zero; \
  939. mkl_dc_type a00 = dc_access(a, lda, 0, 0); \
  940. mkl_dc_type a10 = dc_access(a, lda, 1, 0); \
  941. mkl_dc_type a20 = dc_access(a, lda, 2, 0); \
  942. mkl_dc_type a01 = dc_access(a, lda, 0, 1); \
  943. mkl_dc_type a11 = dc_access(a, lda, 1, 1); \
  944. mkl_dc_type a21 = dc_access(a, lda, 2, 1); \
  945. mkl_dc_type a02 = dc_access(a, lda, 0, 2); \
  946. mkl_dc_type a12 = dc_access(a, lda, 1, 2); \
  947. mkl_dc_type a22 = dc_access(a, lda, 2, 2); \
  948. MKL_DC_IS_ZERO2(is_zero, a00); \
  949. if (!is_zero) { \
  950. MKL_DC_DIV(a10, a10, a00); \
  951. MKL_DC_DIV(a20, a20, a00); \
  952. } else { \
  953. *info = 1; \
  954. } \
  955. MKL_DC_SUB_MUL(a11, a11, a10, a01); \
  956. MKL_DC_SUB_MUL(a21, a21, a20, a01); \
  957. MKL_DC_SUB_MUL(a12, a12, a10, a02); \
  958. MKL_DC_SUB_MUL(a22, a22, a20, a02); \
  959. MKL_DC_IS_ZERO2(is_zero, a11); \
  960. if (!is_zero) { \
  961. MKL_DC_DIV(a21, a21, a11); \
  962. } else if (*info == 0) { \
  963. *info = 2; \
  964. } \
  965. MKL_DC_SUB_MUL(a22, a22, a21, a12); \
  966. MKL_DC_IS_ZERO2(is_zero, a22); \
  967. if (is_zero && *info == 0) { \
  968. *info = 3; \
  969. } \
  970. dc_access(a, lda, 1, 0) = a10; \
  971. dc_access(a, lda, 2, 0) = a20; \
  972. dc_access(a, lda, 1, 1) = a11; \
  973. dc_access(a, lda, 2, 1) = a21; \
  974. dc_access(a, lda, 1, 2) = a12; \
  975. dc_access(a, lda, 2, 2) = a22; \
  976. } while (0) \
  977. #define mkl_dc_getrfnp_4x4(m, n, a, lda, info, dc_access) \
  978. do { \
  979. MKL_INT is_zero; \
  980. mkl_dc_type a00 = dc_access(a, lda, 0, 0); \
  981. mkl_dc_type a10 = dc_access(a, lda, 1, 0); \
  982. mkl_dc_type a20 = dc_access(a, lda, 2, 0); \
  983. mkl_dc_type a30 = dc_access(a, lda, 3, 0); \
  984. mkl_dc_type a01 = dc_access(a, lda, 0, 1); \
  985. mkl_dc_type a11 = dc_access(a, lda, 1, 1); \
  986. mkl_dc_type a21 = dc_access(a, lda, 2, 1); \
  987. mkl_dc_type a31 = dc_access(a, lda, 3, 1); \
  988. mkl_dc_type a02 = dc_access(a, lda, 0, 2); \
  989. mkl_dc_type a12 = dc_access(a, lda, 1, 2); \
  990. mkl_dc_type a22 = dc_access(a, lda, 2, 2); \
  991. mkl_dc_type a32 = dc_access(a, lda, 3, 2); \
  992. mkl_dc_type a03 = dc_access(a, lda, 0, 3); \
  993. mkl_dc_type a13 = dc_access(a, lda, 1, 3); \
  994. mkl_dc_type a23 = dc_access(a, lda, 2, 3); \
  995. mkl_dc_type a33 = dc_access(a, lda, 3, 3); \
  996. MKL_DC_IS_ZERO2(is_zero, a00); \
  997. if (!is_zero) { \
  998. MKL_DC_DIV(a10, a10, a00); \
  999. MKL_DC_DIV(a20, a20, a00); \
  1000. MKL_DC_DIV(a30, a30, a00); \
  1001. } else { \
  1002. *info = 1; \
  1003. } \
  1004. MKL_DC_SUB_MUL(a11, a11, a10, a01); \
  1005. MKL_DC_SUB_MUL(a21, a21, a20, a01); \
  1006. MKL_DC_SUB_MUL(a31, a31, a30, a01); \
  1007. MKL_DC_SUB_MUL(a12, a12, a10, a02); \
  1008. MKL_DC_SUB_MUL(a22, a22, a20, a02); \
  1009. MKL_DC_SUB_MUL(a32, a32, a30, a02); \
  1010. MKL_DC_SUB_MUL(a13, a13, a10, a03); \
  1011. MKL_DC_SUB_MUL(a23, a23, a20, a03); \
  1012. MKL_DC_SUB_MUL(a33, a33, a30, a03); \
  1013. MKL_DC_IS_ZERO2(is_zero, a11); \
  1014. if (!is_zero) { \
  1015. MKL_DC_DIV(a21, a21, a11); \
  1016. MKL_DC_DIV(a31, a31, a11); \
  1017. } else if (*info == 0) { \
  1018. *info = 2; \
  1019. } \
  1020. MKL_DC_SUB_MUL(a22, a22, a21, a12); \
  1021. MKL_DC_SUB_MUL(a32, a32, a31, a12); \
  1022. MKL_DC_SUB_MUL(a23, a23, a21, a13); \
  1023. MKL_DC_SUB_MUL(a33, a33, a31, a13); \
  1024. MKL_DC_IS_ZERO2(is_zero, a22); \
  1025. if (!is_zero) { \
  1026. MKL_DC_DIV(a32, a32, a22); \
  1027. } else if (*info == 0) { \
  1028. *info = 3; \
  1029. } \
  1030. MKL_DC_SUB_MUL(a33, a33, a32, a23); \
  1031. MKL_DC_IS_ZERO2(is_zero, a33); \
  1032. if (is_zero && *info == 0) { \
  1033. *info = 4; \
  1034. } \
  1035. dc_access(a, lda, 1, 0) = a10; \
  1036. dc_access(a, lda, 2, 0) = a20; \
  1037. dc_access(a, lda, 3, 0) = a30; \
  1038. dc_access(a, lda, 1, 1) = a11; \
  1039. dc_access(a, lda, 2, 1) = a21; \
  1040. dc_access(a, lda, 3, 1) = a31; \
  1041. dc_access(a, lda, 1, 2) = a12; \
  1042. dc_access(a, lda, 2, 2) = a22; \
  1043. dc_access(a, lda, 3, 2) = a32; \
  1044. dc_access(a, lda, 1, 3) = a13; \
  1045. dc_access(a, lda, 2, 3) = a23; \
  1046. dc_access(a, lda, 3, 3) = a33; \
  1047. } while (0)
  1048. #define mkl_dc_getrfnp_5x5(m, n, a, lda, info, dc_access) \
  1049. do { \
  1050. MKL_INT is_zero; \
  1051. mkl_dc_type a00 = dc_access(a, lda, 0, 0); \
  1052. mkl_dc_type a10 = dc_access(a, lda, 1, 0); \
  1053. mkl_dc_type a20 = dc_access(a, lda, 2, 0); \
  1054. mkl_dc_type a30 = dc_access(a, lda, 3, 0); \
  1055. mkl_dc_type a40 = dc_access(a, lda, 4, 0); \
  1056. mkl_dc_type a01 = dc_access(a, lda, 0, 1); \
  1057. mkl_dc_type a11 = dc_access(a, lda, 1, 1); \
  1058. mkl_dc_type a21 = dc_access(a, lda, 2, 1); \
  1059. mkl_dc_type a31 = dc_access(a, lda, 3, 1); \
  1060. mkl_dc_type a41 = dc_access(a, lda, 4, 1); \
  1061. mkl_dc_type a02 = dc_access(a, lda, 0, 2); \
  1062. mkl_dc_type a12 = dc_access(a, lda, 1, 2); \
  1063. mkl_dc_type a22 = dc_access(a, lda, 2, 2); \
  1064. mkl_dc_type a32 = dc_access(a, lda, 3, 2); \
  1065. mkl_dc_type a42 = dc_access(a, lda, 4, 2); \
  1066. mkl_dc_type a03 = dc_access(a, lda, 0, 3); \
  1067. mkl_dc_type a13 = dc_access(a, lda, 1, 3); \
  1068. mkl_dc_type a23 = dc_access(a, lda, 2, 3); \
  1069. mkl_dc_type a33 = dc_access(a, lda, 3, 3); \
  1070. mkl_dc_type a43 = dc_access(a, lda, 4, 3); \
  1071. mkl_dc_type a04 = dc_access(a, lda, 0, 4); \
  1072. mkl_dc_type a14 = dc_access(a, lda, 1, 4); \
  1073. mkl_dc_type a24 = dc_access(a, lda, 2, 4); \
  1074. mkl_dc_type a34 = dc_access(a, lda, 3, 4); \
  1075. mkl_dc_type a44 = dc_access(a, lda, 4, 4); \
  1076. MKL_DC_IS_ZERO2(is_zero, a00); \
  1077. if (!is_zero) { \
  1078. MKL_DC_DIV(a10, a10, a00); \
  1079. MKL_DC_DIV(a20, a20, a00); \
  1080. MKL_DC_DIV(a30, a30, a00); \
  1081. MKL_DC_DIV(a40, a40, a00); \
  1082. } else { \
  1083. *info = 1; \
  1084. } \
  1085. MKL_DC_SUB_MUL(a11, a11, a10, a01); \
  1086. MKL_DC_SUB_MUL(a21, a21, a20, a01); \
  1087. MKL_DC_SUB_MUL(a31, a31, a30, a01); \
  1088. MKL_DC_SUB_MUL(a41, a41, a40, a01); \
  1089. MKL_DC_SUB_MUL(a12, a12, a10, a02); \
  1090. MKL_DC_SUB_MUL(a22, a22, a20, a02); \
  1091. MKL_DC_SUB_MUL(a32, a32, a30, a02); \
  1092. MKL_DC_SUB_MUL(a42, a42, a40, a02); \
  1093. MKL_DC_SUB_MUL(a13, a13, a10, a03); \
  1094. MKL_DC_SUB_MUL(a23, a23, a20, a03); \
  1095. MKL_DC_SUB_MUL(a33, a33, a30, a03); \
  1096. MKL_DC_SUB_MUL(a43, a43, a40, a03); \
  1097. MKL_DC_SUB_MUL(a14, a14, a10, a04); \
  1098. MKL_DC_SUB_MUL(a24, a24, a20, a04); \
  1099. MKL_DC_SUB_MUL(a34, a34, a30, a04); \
  1100. MKL_DC_SUB_MUL(a44, a44, a40, a04); \
  1101. MKL_DC_IS_ZERO2(is_zero, a11); \
  1102. if (!is_zero) { \
  1103. MKL_DC_DIV(a21, a21, a11); \
  1104. MKL_DC_DIV(a31, a31, a11); \
  1105. MKL_DC_DIV(a41, a41, a11); \
  1106. } else if (*info == 0) { \
  1107. *info = 2; \
  1108. } \
  1109. MKL_DC_SUB_MUL(a22, a22, a21, a12); \
  1110. MKL_DC_SUB_MUL(a32, a32, a31, a12); \
  1111. MKL_DC_SUB_MUL(a42, a42, a41, a12); \
  1112. MKL_DC_SUB_MUL(a23, a23, a21, a13); \
  1113. MKL_DC_SUB_MUL(a33, a33, a31, a13); \
  1114. MKL_DC_SUB_MUL(a43, a43, a41, a13); \
  1115. MKL_DC_SUB_MUL(a24, a24, a21, a14); \
  1116. MKL_DC_SUB_MUL(a34, a34, a31, a14); \
  1117. MKL_DC_SUB_MUL(a44, a44, a41, a14); \
  1118. MKL_DC_IS_ZERO2(is_zero, a22); \
  1119. if (!is_zero) { \
  1120. MKL_DC_DIV(a32, a32, a22); \
  1121. MKL_DC_DIV(a42, a42, a22); \
  1122. } else if (*info == 0) { \
  1123. *info = 3; \
  1124. } \
  1125. MKL_DC_SUB_MUL(a33, a33, a32, a23); \
  1126. MKL_DC_SUB_MUL(a43, a43, a42, a23); \
  1127. MKL_DC_SUB_MUL(a34, a34, a32, a24); \
  1128. MKL_DC_SUB_MUL(a44, a44, a42, a24); \
  1129. MKL_DC_IS_ZERO2(is_zero, a33); \
  1130. if (!is_zero) { \
  1131. MKL_DC_DIV(a43, a43, a33); \
  1132. } else if (*info == 0) { \
  1133. *info = 4; \
  1134. } \
  1135. MKL_DC_SUB_MUL(a44, a44, a43, a34); \
  1136. MKL_DC_IS_ZERO2(is_zero, a44); \
  1137. if (is_zero && *info == 0) { \
  1138. *info = 5; \
  1139. } \
  1140. dc_access(a, lda, 1, 0) = a10; \
  1141. dc_access(a, lda, 2, 0) = a20; \
  1142. dc_access(a, lda, 3, 0) = a30; \
  1143. dc_access(a, lda, 4, 0) = a40; \
  1144. dc_access(a, lda, 1, 1) = a11; \
  1145. dc_access(a, lda, 2, 1) = a21; \
  1146. dc_access(a, lda, 3, 1) = a31; \
  1147. dc_access(a, lda, 4, 1) = a41; \
  1148. dc_access(a, lda, 1, 2) = a12; \
  1149. dc_access(a, lda, 2, 2) = a22; \
  1150. dc_access(a, lda, 3, 2) = a32; \
  1151. dc_access(a, lda, 4, 2) = a42; \
  1152. dc_access(a, lda, 1, 3) = a13; \
  1153. dc_access(a, lda, 2, 3) = a23; \
  1154. dc_access(a, lda, 3, 3) = a33; \
  1155. dc_access(a, lda, 4, 3) = a43; \
  1156. dc_access(a, lda, 1, 4) = a14; \
  1157. dc_access(a, lda, 2, 4) = a24; \
  1158. dc_access(a, lda, 3, 4) = a34; \
  1159. dc_access(a, lda, 4, 4) = a44; \
  1160. } while (0)
  1161. #define mkl_dc_getrfnp_generic(m, n, a, lda, info, dc_access) \
  1162. do { \
  1163. if (m == n && m <= 5) { \
  1164. switch (m) { \
  1165. case 0: \
  1166. break; \
  1167. case 1: \
  1168. mkl_dc_getrfnp_1x1(m, n, a, lda, info, dc_access); \
  1169. break; \
  1170. case 2: \
  1171. mkl_dc_getrfnp_2x2(m, n, a, lda, info, dc_access); \
  1172. break; \
  1173. case 3: \
  1174. mkl_dc_getrfnp_3x3(m, n, a, lda, info, dc_access); \
  1175. break; \
  1176. case 4: \
  1177. mkl_dc_getrfnp_4x4(m, n, a, lda, info, dc_access); \
  1178. break; \
  1179. case 5: \
  1180. mkl_dc_getrfnp_5x5(m, n, a, lda, info, dc_access); \
  1181. break; \
  1182. } \
  1183. } else { \
  1184. mkl_dc_getrfnp_mxn(m, n, a, lda, info, dc_access); \
  1185. } \
  1186. } while (0)
  1187. static __inline void mkl_dc_getrfnp(MKL_INT m, MKL_INT n, mkl_dc_type* a, MKL_INT lda, MKL_INT* info)
  1188. {
  1189. *info = 0;
  1190. mkl_dc_getrfnp_generic(m, n, a, lda, info, MKL_DC_MN);
  1191. }
  1192. #ifndef MKL_DIRECT_CALL_LAPACKE_DISABLE
  1193. static __inline lapack_int mkl_dc_lapacke_getrfnp_convert(int matrix_layout, lapack_int m, lapack_int n, mkl_dc_type* a, lapack_int lda)
  1194. {
  1195. lapack_int info = 0;
  1196. if (MKL_DC_GETRFNP_CHECKSIZE(m, n)) {
  1197. if (matrix_layout == LAPACK_ROW_MAJOR) {
  1198. mkl_dc_getrfnp_generic(m, n, a, lda, &info, MKL_DC_MT);
  1199. } else {
  1200. mkl_dc_getrfnp_generic(m, n, a, lda, &info, MKL_DC_MN);
  1201. }
  1202. return info;
  1203. }
  1204. return MKL_DC_CONCAT3(LAPACKE_mkl_, MKL_DC_PREC_LETTER, getrfnp)(matrix_layout, m, n, a, lda);
  1205. }
  1206. #endif
  1207. /* ?GETRS */
  1208. #define mkl_dc_getrs_nxnrhs(trans, n, nrhs, a, lda, ipiv, b, ldb, info, dc_access, b_access) \
  1209. do { \
  1210. MKL_INT i, j, k; \
  1211. if (trans == 'N') { \
  1212. for (i = 0; i < n; ++i) { \
  1213. MKL_INT ip = ipiv[i] - 1; \
  1214. if (ip != i) { \
  1215. for (k = 0; k < nrhs; k++) { \
  1216. MKL_DC_SWAP(b_access(b, ldb, i, k), b_access(b, ldb, ip, k)); \
  1217. } \
  1218. } \
  1219. } \
  1220. for (j = 0; j < nrhs; ++j) { \
  1221. for (k = 0; k < n; k++) { \
  1222. if (!MKL_DC_IS_ZERO(b_access(b, ldb, k, j))) { \
  1223. for (i = k + 1; i < n; ++i) { \
  1224. MKL_DC_SUB_MUL(b_access(b, ldb, i, j), b_access(b, ldb, i, j), b_access(b, ldb, k, j), dc_access(a, lda, i, k)); \
  1225. } \
  1226. } \
  1227. } \
  1228. } \
  1229. for (j = 0; j < nrhs; ++j) { \
  1230. for (k = n - 1; k >= 0; k--) { \
  1231. if (!MKL_DC_IS_ZERO(b_access(b, ldb, k, j))) { \
  1232. MKL_DC_DIV(b_access(b, ldb, k, j), b_access(b, ldb, k, j), dc_access(a, lda, k, k)); \
  1233. for (i = 0; i < k; ++i) { \
  1234. MKL_DC_SUB_MUL(b_access(b, ldb, i, j), b_access(b, ldb, i, j), b_access(b, ldb, k, j), dc_access(a, lda, i, k)); \
  1235. } \
  1236. } \
  1237. } \
  1238. } \
  1239. } else { \
  1240. for (j = 0; j < nrhs; ++j) { \
  1241. for (i = 0; i < n; ++i) { \
  1242. mkl_dc_type temp = b_access(b, ldb, i, j); \
  1243. mkl_dc_type aii = dc_access(a, lda, i, i); \
  1244. for (k = 0; k < i; k++) { \
  1245. mkl_dc_type aki = dc_access(a, lda, k, i); \
  1246. if (trans == 'C') { \
  1247. MKL_DC_CONJ(aki, aki); \
  1248. } \
  1249. MKL_DC_SUB_MUL(temp, temp, aki, b_access(b, ldb, k, j)); \
  1250. } \
  1251. if (trans == 'C') { \
  1252. MKL_DC_CONJ(aii, aii); \
  1253. } \
  1254. MKL_DC_DIV(b_access(b, ldb, i, j), temp, aii); \
  1255. } \
  1256. } \
  1257. for (j = 0; j < nrhs; ++j) { \
  1258. for (i = n - 1; i >= 0; i--) { \
  1259. mkl_dc_type temp = b_access(b, ldb, i, j); \
  1260. for (k = i + 1; k < n; k++) { \
  1261. mkl_dc_type aki = dc_access(a, lda, k, i); \
  1262. if (trans == 'C') { \
  1263. MKL_DC_CONJ(aki, aki); \
  1264. } \
  1265. MKL_DC_SUB_MUL(temp, temp, aki, b_access(b, ldb, k, j)); \
  1266. } \
  1267. b_access(b, ldb, i, j) = temp; \
  1268. } \
  1269. } \
  1270. for (i = n - 1; i >= 0; i--) { \
  1271. MKL_INT ip = ipiv[i] - 1; \
  1272. if (ip != i) { \
  1273. for (k = 0; k < nrhs; k++) { \
  1274. MKL_DC_SWAP(b_access(b, ldb, i, k), b_access(b, ldb, ip, k)); \
  1275. } \
  1276. } \
  1277. } \
  1278. } \
  1279. } while (0)
  1280. #define mkl_dc_getrs_trans(trans, n, nrhs, a, lda, ipiv, b, ldb, info, dc_access, b_access) \
  1281. do { \
  1282. switch (n) { \
  1283. case 1: \
  1284. mkl_dc_getrs_nxnrhs(trans, 1, nrhs, a, lda, ipiv, b, ldb, info, dc_access, b_access); \
  1285. break; \
  1286. case 2: \
  1287. mkl_dc_getrs_nxnrhs(trans, 2, nrhs, a, lda, ipiv, b, ldb, info, dc_access, b_access); \
  1288. break; \
  1289. case 3: \
  1290. mkl_dc_getrs_nxnrhs(trans, 3, nrhs, a, lda, ipiv, b, ldb, info, dc_access, b_access); \
  1291. break; \
  1292. case 4: \
  1293. mkl_dc_getrs_nxnrhs(trans, 4, nrhs, a, lda, ipiv, b, ldb, info, dc_access, b_access); \
  1294. break; \
  1295. case 5: \
  1296. mkl_dc_getrs_nxnrhs(trans, 5, nrhs, a, lda, ipiv, b, ldb, info, dc_access, b_access); \
  1297. break; \
  1298. default: \
  1299. mkl_dc_getrs_nxnrhs(trans, n, nrhs, a, lda, ipiv, b, ldb, info, dc_access, b_access); \
  1300. break; \
  1301. } \
  1302. } while (0)
  1303. #define mkl_dc_getrs_generic(trans, n, nrhs, a, lda, ipiv, b, ldb, info, dc_access, b_access) \
  1304. do { \
  1305. int ntrans = MKL_DC_MisN(trans); \
  1306. int ttrans = MKL_DC_MisT(trans); \
  1307. *info = 0; \
  1308. if (n != 0 && nrhs != 0) { \
  1309. if (ntrans) { \
  1310. mkl_dc_getrs_trans('N', n, nrhs, a, lda, ipiv, b, ldb, info, dc_access, b_access); \
  1311. } else if (ttrans) { \
  1312. mkl_dc_getrs_trans('T', n, nrhs, a, lda, ipiv, b, ldb, info, dc_access, b_access); \
  1313. } else { \
  1314. mkl_dc_getrs_trans('C', n, nrhs, a, lda, ipiv, b, ldb, info, dc_access, b_access); \
  1315. } \
  1316. } \
  1317. } while (0)
  1318. static __inline void mkl_dc_getrs(char trans, MKL_INT n, MKL_INT nrhs, const mkl_dc_type* a, MKL_INT lda,
  1319. const MKL_INT* ipiv, mkl_dc_type* b, MKL_INT ldb, MKL_INT* info)
  1320. {
  1321. mkl_dc_getrs_generic(trans, n, nrhs, a, lda, ipiv, b, ldb, info, MKL_DC_MN, MKL_DC_MN);
  1322. }
  1323. #ifndef MKL_DIRECT_CALL_LAPACKE_DISABLE
  1324. static __inline lapack_int mkl_dc_lapacke_getrs_convert(int matrix_layout, char trans, lapack_int n, lapack_int nrhs,
  1325. const mkl_dc_type* a, lapack_int lda, const lapack_int* ipiv,
  1326. mkl_dc_type* b, lapack_int ldb)
  1327. {
  1328. lapack_int info = 0;
  1329. if (MKL_DC_GETRS_CHECKSIZE(n, nrhs)) {
  1330. if (matrix_layout == LAPACK_ROW_MAJOR) {
  1331. mkl_dc_getrs_generic(trans, n, nrhs, a, lda, ipiv, b, ldb, &info, MKL_DC_MT, MKL_DC_MT);
  1332. } else {
  1333. mkl_dc_getrs_generic(trans, n, nrhs, a, lda, ipiv, b, ldb, &info, MKL_DC_MN, MKL_DC_MN);
  1334. }
  1335. return info;
  1336. }
  1337. return MKL_DC_CONCAT3(LAPACKE_, MKL_DC_PREC_LETTER, getrs)(matrix_layout, trans, n, nrhs, a, lda, ipiv, b, ldb);
  1338. }
  1339. #endif
  1340. /* ?GETRI */
  1341. #define mkl_dc_getri_1x1(n, a, lda, ipiv, info, dc_access) \
  1342. do { \
  1343. mkl_dc_type a00 = dc_access(a, lda, 0, 0); \
  1344. mkl_dc_type ajj; \
  1345. if (MKL_DC_IS_ZERO(a00)) { \
  1346. *info = 1; \
  1347. break; \
  1348. } \
  1349. *info = 0; \
  1350. MKL_DC_INV(a00, a00); \
  1351. MKL_DC_NEG(ajj, a00); \
  1352. dc_access(a, lda, 0, 0) = a00; \
  1353. } while (0)
  1354. #define mkl_dc_getri_2x2(n, a, lda, ipiv, info, dc_access) \
  1355. do { \
  1356. mkl_dc_type a00 = dc_access(a, lda, 0, 0); \
  1357. mkl_dc_type a01 = dc_access(a, lda, 0, 1); \
  1358. mkl_dc_type a10 = dc_access(a, lda, 1, 0); \
  1359. mkl_dc_type a11 = dc_access(a, lda, 1, 1); \
  1360. mkl_dc_type w1; \
  1361. mkl_dc_type tmp; \
  1362. MKL_INT jp; \
  1363. mkl_dc_type ajj; \
  1364. if (MKL_DC_IS_ZERO(a00)) { \
  1365. *info = 1; \
  1366. break; \
  1367. } \
  1368. if (MKL_DC_IS_ZERO(a11)) { \
  1369. *info = 2; \
  1370. break; \
  1371. } \
  1372. *info = 0; \
  1373. MKL_DC_INV(a00, a00); \
  1374. MKL_DC_NEG(ajj, a00); \
  1375. MKL_DC_INV(a11, a11); \
  1376. MKL_DC_NEG(ajj, a11); \
  1377. tmp = a01; \
  1378. MKL_DC_MUL(a01, a01, a00); \
  1379. MKL_DC_MUL(a01, a01, ajj); \
  1380. w1 = a10; \
  1381. MKL_DC_SET_ZERO(a10); \
  1382. MKL_DC_NEG(tmp, w1); \
  1383. MKL_DC_MUL_ADD(a00, tmp, a01, a00); \
  1384. MKL_DC_MUL_ADD(a10, tmp, a11, a10); \
  1385. dc_access(a, lda, 0, 0) = a00; \
  1386. dc_access(a, lda, 0, 1) = a01; \
  1387. dc_access(a, lda, 1, 0) = a10; \
  1388. dc_access(a, lda, 1, 1) = a11; \
  1389. jp = ipiv[0] - 1; \
  1390. if (jp != 0) { \
  1391. MKL_DC_SWAP(dc_access(a, lda, 0, 0), dc_access(a, lda, 0, jp)); \
  1392. MKL_DC_SWAP(dc_access(a, lda, 1, 0), dc_access(a, lda, 1, jp)); \
  1393. } \
  1394. } while (0)
  1395. #define mkl_dc_getri_3x3(n, a, lda, ipiv, info, dc_access) \
  1396. do { \
  1397. mkl_dc_type a00 = dc_access(a, lda, 0, 0); \
  1398. mkl_dc_type a01 = dc_access(a, lda, 0, 1); \
  1399. mkl_dc_type a02 = dc_access(a, lda, 0, 2); \
  1400. mkl_dc_type a10 = dc_access(a, lda, 1, 0); \
  1401. mkl_dc_type a11 = dc_access(a, lda, 1, 1); \
  1402. mkl_dc_type a12 = dc_access(a, lda, 1, 2); \
  1403. mkl_dc_type a20 = dc_access(a, lda, 2, 0); \
  1404. mkl_dc_type a21 = dc_access(a, lda, 2, 1); \
  1405. mkl_dc_type a22 = dc_access(a, lda, 2, 2); \
  1406. mkl_dc_type w1; \
  1407. mkl_dc_type w2; \
  1408. mkl_dc_type tmp; \
  1409. MKL_INT jp; \
  1410. mkl_dc_type ajj; \
  1411. if (MKL_DC_IS_ZERO(a00)) { \
  1412. *info = 1; \
  1413. break; \
  1414. } \
  1415. if (MKL_DC_IS_ZERO(a11)) { \
  1416. *info = 2; \
  1417. break; \
  1418. } \
  1419. if (MKL_DC_IS_ZERO(a22)) { \
  1420. *info = 3; \
  1421. break; \
  1422. } \
  1423. *info = 0; \
  1424. MKL_DC_INV(a00, a00); \
  1425. MKL_DC_NEG(ajj, a00); \
  1426. MKL_DC_INV(a11, a11); \
  1427. MKL_DC_NEG(ajj, a11); \
  1428. tmp = a01; \
  1429. MKL_DC_MUL(a01, a01, a00); \
  1430. MKL_DC_MUL(a01, a01, ajj); \
  1431. MKL_DC_INV(a22, a22); \
  1432. MKL_DC_NEG(ajj, a22); \
  1433. tmp = a02; \
  1434. MKL_DC_MUL(a02, a02, a00); \
  1435. tmp = a12; \
  1436. MKL_DC_MUL_ADD(a02, tmp, a01, a02); \
  1437. MKL_DC_MUL(a12, a12, a11); \
  1438. MKL_DC_MUL(a02, a02, ajj); \
  1439. MKL_DC_MUL(a12, a12, ajj); \
  1440. w2 = a21; \
  1441. MKL_DC_SET_ZERO(a21); \
  1442. MKL_DC_NEG(tmp, w2); \
  1443. MKL_DC_MUL_ADD(a01, tmp, a02, a01); \
  1444. MKL_DC_MUL_ADD(a11, tmp, a12, a11); \
  1445. MKL_DC_MUL_ADD(a21, tmp, a22, a21); \
  1446. w1 = a10; \
  1447. MKL_DC_SET_ZERO(a10); \
  1448. w2 = a20; \
  1449. MKL_DC_SET_ZERO(a20); \
  1450. MKL_DC_NEG(tmp, w1); \
  1451. MKL_DC_MUL_ADD(a00, tmp, a01, a00); \
  1452. MKL_DC_MUL_ADD(a10, tmp, a11, a10); \
  1453. MKL_DC_MUL_ADD(a20, tmp, a21, a20); \
  1454. MKL_DC_NEG(tmp, w2); \
  1455. MKL_DC_MUL_ADD(a00, tmp, a02, a00); \
  1456. MKL_DC_MUL_ADD(a10, tmp, a12, a10); \
  1457. MKL_DC_MUL_ADD(a20, tmp, a22, a20); \
  1458. dc_access(a, lda, 0, 0) = a00; \
  1459. dc_access(a, lda, 0, 1) = a01; \
  1460. dc_access(a, lda, 0, 2) = a02; \
  1461. dc_access(a, lda, 1, 0) = a10; \
  1462. dc_access(a, lda, 1, 1) = a11; \
  1463. dc_access(a, lda, 1, 2) = a12; \
  1464. dc_access(a, lda, 2, 0) = a20; \
  1465. dc_access(a, lda, 2, 1) = a21; \
  1466. dc_access(a, lda, 2, 2) = a22; \
  1467. jp = ipiv[1] - 1; \
  1468. if (jp != 1) { \
  1469. MKL_DC_SWAP(dc_access(a, lda, 0, 1), dc_access(a, lda, 0, jp)); \
  1470. MKL_DC_SWAP(dc_access(a, lda, 1, 1), dc_access(a, lda, 1, jp)); \
  1471. MKL_DC_SWAP(dc_access(a, lda, 2, 1), dc_access(a, lda, 2, jp)); \
  1472. } \
  1473. jp = ipiv[0] - 1; \
  1474. if (jp != 0) { \
  1475. MKL_DC_SWAP(dc_access(a, lda, 0, 0), dc_access(a, lda, 0, jp)); \
  1476. MKL_DC_SWAP(dc_access(a, lda, 1, 0), dc_access(a, lda, 1, jp)); \
  1477. MKL_DC_SWAP(dc_access(a, lda, 2, 0), dc_access(a, lda, 2, jp)); \
  1478. } \
  1479. } while (0)
  1480. #define mkl_dc_getri_4x4(n, a, lda, ipiv, info, dc_access) \
  1481. do { \
  1482. mkl_dc_type a00 = dc_access(a, lda, 0, 0); \
  1483. mkl_dc_type a01 = dc_access(a, lda, 0, 1); \
  1484. mkl_dc_type a02 = dc_access(a, lda, 0, 2); \
  1485. mkl_dc_type a03 = dc_access(a, lda, 0, 3); \
  1486. mkl_dc_type a10 = dc_access(a, lda, 1, 0); \
  1487. mkl_dc_type a11 = dc_access(a, lda, 1, 1); \
  1488. mkl_dc_type a12 = dc_access(a, lda, 1, 2); \
  1489. mkl_dc_type a13 = dc_access(a, lda, 1, 3); \
  1490. mkl_dc_type a20 = dc_access(a, lda, 2, 0); \
  1491. mkl_dc_type a21 = dc_access(a, lda, 2, 1); \
  1492. mkl_dc_type a22 = dc_access(a, lda, 2, 2); \
  1493. mkl_dc_type a23 = dc_access(a, lda, 2, 3); \
  1494. mkl_dc_type a30 = dc_access(a, lda, 3, 0); \
  1495. mkl_dc_type a31 = dc_access(a, lda, 3, 1); \
  1496. mkl_dc_type a32 = dc_access(a, lda, 3, 2); \
  1497. mkl_dc_type a33 = dc_access(a, lda, 3, 3); \
  1498. mkl_dc_type w1; \
  1499. mkl_dc_type w2; \
  1500. mkl_dc_type w3; \
  1501. mkl_dc_type tmp; \
  1502. MKL_INT jp; \
  1503. mkl_dc_type ajj; \
  1504. if (MKL_DC_IS_ZERO(a00)) { \
  1505. *info = 1; \
  1506. break; \
  1507. } \
  1508. if (MKL_DC_IS_ZERO(a11)) { \
  1509. *info = 2; \
  1510. break; \
  1511. } \
  1512. if (MKL_DC_IS_ZERO(a22)) { \
  1513. *info = 3; \
  1514. break; \
  1515. } \
  1516. if (MKL_DC_IS_ZERO(a33)) { \
  1517. *info = 4; \
  1518. break; \
  1519. } \
  1520. *info = 0; \
  1521. MKL_DC_INV(a00, a00); \
  1522. MKL_DC_NEG(ajj, a00); \
  1523. MKL_DC_INV(a11, a11); \
  1524. MKL_DC_NEG(ajj, a11); \
  1525. tmp = a01; \
  1526. MKL_DC_MUL(a01, a01, a00); \
  1527. MKL_DC_MUL(a01, a01, ajj); \
  1528. MKL_DC_INV(a22, a22); \
  1529. MKL_DC_NEG(ajj, a22); \
  1530. tmp = a02; \
  1531. MKL_DC_MUL(a02, a02, a00); \
  1532. tmp = a12; \
  1533. MKL_DC_MUL_ADD(a02, tmp, a01, a02); \
  1534. MKL_DC_MUL(a12, a12, a11); \
  1535. MKL_DC_MUL(a02, a02, ajj); \
  1536. MKL_DC_MUL(a12, a12, ajj); \
  1537. MKL_DC_INV(a33, a33); \
  1538. MKL_DC_NEG(ajj, a33); \
  1539. tmp = a03; \
  1540. MKL_DC_MUL(a03, a03, a00); \
  1541. tmp = a13; \
  1542. MKL_DC_MUL_ADD(a03, tmp, a01, a03); \
  1543. MKL_DC_MUL(a13, a13, a11); \
  1544. tmp = a23; \
  1545. MKL_DC_MUL_ADD(a03, tmp, a02, a03); \
  1546. MKL_DC_MUL_ADD(a13, tmp, a12, a13); \
  1547. MKL_DC_MUL(a23, a23, a22); \
  1548. MKL_DC_MUL(a03, a03, ajj); \
  1549. MKL_DC_MUL(a13, a13, ajj); \
  1550. MKL_DC_MUL(a23, a23, ajj); \
  1551. w3 = a32; \
  1552. MKL_DC_SET_ZERO(a32); \
  1553. MKL_DC_NEG(tmp, w3); \
  1554. MKL_DC_MUL_ADD(a02, tmp, a03, a02); \
  1555. MKL_DC_MUL_ADD(a12, tmp, a13, a12); \
  1556. MKL_DC_MUL_ADD(a22, tmp, a23, a22); \
  1557. MKL_DC_MUL_ADD(a32, tmp, a33, a32); \
  1558. w2 = a21; \
  1559. MKL_DC_SET_ZERO(a21); \
  1560. w3 = a31; \
  1561. MKL_DC_SET_ZERO(a31); \
  1562. MKL_DC_NEG(tmp, w2); \
  1563. MKL_DC_MUL_ADD(a01, tmp, a02, a01); \
  1564. MKL_DC_MUL_ADD(a11, tmp, a12, a11); \
  1565. MKL_DC_MUL_ADD(a21, tmp, a22, a21); \
  1566. MKL_DC_MUL_ADD(a31, tmp, a32, a31); \
  1567. MKL_DC_NEG(tmp, w3); \
  1568. MKL_DC_MUL_ADD(a01, tmp, a03, a01); \
  1569. MKL_DC_MUL_ADD(a11, tmp, a13, a11); \
  1570. MKL_DC_MUL_ADD(a21, tmp, a23, a21); \
  1571. MKL_DC_MUL_ADD(a31, tmp, a33, a31); \
  1572. w1 = a10; \
  1573. MKL_DC_SET_ZERO(a10); \
  1574. w2 = a20; \
  1575. MKL_DC_SET_ZERO(a20); \
  1576. w3 = a30; \
  1577. MKL_DC_SET_ZERO(a30); \
  1578. MKL_DC_NEG(tmp, w1); \
  1579. MKL_DC_MUL_ADD(a00, tmp, a01, a00); \
  1580. MKL_DC_MUL_ADD(a10, tmp, a11, a10); \
  1581. MKL_DC_MUL_ADD(a20, tmp, a21, a20); \
  1582. MKL_DC_MUL_ADD(a30, tmp, a31, a30); \
  1583. MKL_DC_NEG(tmp, w2); \
  1584. MKL_DC_MUL_ADD(a00, tmp, a02, a00); \
  1585. MKL_DC_MUL_ADD(a10, tmp, a12, a10); \
  1586. MKL_DC_MUL_ADD(a20, tmp, a22, a20); \
  1587. MKL_DC_MUL_ADD(a30, tmp, a32, a30); \
  1588. MKL_DC_NEG(tmp, w3); \
  1589. MKL_DC_MUL_ADD(a00, tmp, a03, a00); \
  1590. MKL_DC_MUL_ADD(a10, tmp, a13, a10); \
  1591. MKL_DC_MUL_ADD(a20, tmp, a23, a20); \
  1592. MKL_DC_MUL_ADD(a30, tmp, a33, a30); \
  1593. dc_access(a, lda, 0, 0) = a00; \
  1594. dc_access(a, lda, 0, 1) = a01; \
  1595. dc_access(a, lda, 0, 2) = a02; \
  1596. dc_access(a, lda, 0, 3) = a03; \
  1597. dc_access(a, lda, 1, 0) = a10; \
  1598. dc_access(a, lda, 1, 1) = a11; \
  1599. dc_access(a, lda, 1, 2) = a12; \
  1600. dc_access(a, lda, 1, 3) = a13; \
  1601. dc_access(a, lda, 2, 0) = a20; \
  1602. dc_access(a, lda, 2, 1) = a21; \
  1603. dc_access(a, lda, 2, 2) = a22; \
  1604. dc_access(a, lda, 2, 3) = a23; \
  1605. dc_access(a, lda, 3, 0) = a30; \
  1606. dc_access(a, lda, 3, 1) = a31; \
  1607. dc_access(a, lda, 3, 2) = a32; \
  1608. dc_access(a, lda, 3, 3) = a33; \
  1609. jp = ipiv[2] - 1; \
  1610. if (jp != 2) { \
  1611. MKL_DC_SWAP(dc_access(a, lda, 0, 2), dc_access(a, lda, 0, jp)); \
  1612. MKL_DC_SWAP(dc_access(a, lda, 1, 2), dc_access(a, lda, 1, jp)); \
  1613. MKL_DC_SWAP(dc_access(a, lda, 2, 2), dc_access(a, lda, 2, jp)); \
  1614. MKL_DC_SWAP(dc_access(a, lda, 3, 2), dc_access(a, lda, 3, jp)); \
  1615. } \
  1616. jp = ipiv[1] - 1; \
  1617. if (jp != 1) { \
  1618. MKL_DC_SWAP(dc_access(a, lda, 0, 1), dc_access(a, lda, 0, jp)); \
  1619. MKL_DC_SWAP(dc_access(a, lda, 1, 1), dc_access(a, lda, 1, jp)); \
  1620. MKL_DC_SWAP(dc_access(a, lda, 2, 1), dc_access(a, lda, 2, jp)); \
  1621. MKL_DC_SWAP(dc_access(a, lda, 3, 1), dc_access(a, lda, 3, jp)); \
  1622. } \
  1623. jp = ipiv[0] - 1; \
  1624. if (jp != 0) { \
  1625. MKL_DC_SWAP(dc_access(a, lda, 0, 0), dc_access(a, lda, 0, jp)); \
  1626. MKL_DC_SWAP(dc_access(a, lda, 1, 0), dc_access(a, lda, 1, jp)); \
  1627. MKL_DC_SWAP(dc_access(a, lda, 2, 0), dc_access(a, lda, 2, jp)); \
  1628. MKL_DC_SWAP(dc_access(a, lda, 3, 0), dc_access(a, lda, 3, jp)); \
  1629. } \
  1630. } while (0)
  1631. #define mkl_dc_getri_5x5(n, a, lda, ipiv, info, dc_access) \
  1632. do { \
  1633. mkl_dc_type a00 = dc_access(a, lda, 0, 0); \
  1634. mkl_dc_type a01 = dc_access(a, lda, 0, 1); \
  1635. mkl_dc_type a02 = dc_access(a, lda, 0, 2); \
  1636. mkl_dc_type a03 = dc_access(a, lda, 0, 3); \
  1637. mkl_dc_type a04 = dc_access(a, lda, 0, 4); \
  1638. mkl_dc_type a10 = dc_access(a, lda, 1, 0); \
  1639. mkl_dc_type a11 = dc_access(a, lda, 1, 1); \
  1640. mkl_dc_type a12 = dc_access(a, lda, 1, 2); \
  1641. mkl_dc_type a13 = dc_access(a, lda, 1, 3); \
  1642. mkl_dc_type a14 = dc_access(a, lda, 1, 4); \
  1643. mkl_dc_type a20 = dc_access(a, lda, 2, 0); \
  1644. mkl_dc_type a21 = dc_access(a, lda, 2, 1); \
  1645. mkl_dc_type a22 = dc_access(a, lda, 2, 2); \
  1646. mkl_dc_type a23 = dc_access(a, lda, 2, 3); \
  1647. mkl_dc_type a24 = dc_access(a, lda, 2, 4); \
  1648. mkl_dc_type a30 = dc_access(a, lda, 3, 0); \
  1649. mkl_dc_type a31 = dc_access(a, lda, 3, 1); \
  1650. mkl_dc_type a32 = dc_access(a, lda, 3, 2); \
  1651. mkl_dc_type a33 = dc_access(a, lda, 3, 3); \
  1652. mkl_dc_type a34 = dc_access(a, lda, 3, 4); \
  1653. mkl_dc_type a40 = dc_access(a, lda, 4, 0); \
  1654. mkl_dc_type a41 = dc_access(a, lda, 4, 1); \
  1655. mkl_dc_type a42 = dc_access(a, lda, 4, 2); \
  1656. mkl_dc_type a43 = dc_access(a, lda, 4, 3); \
  1657. mkl_dc_type a44 = dc_access(a, lda, 4, 4); \
  1658. mkl_dc_type w1; \
  1659. mkl_dc_type w2; \
  1660. mkl_dc_type w3; \
  1661. mkl_dc_type w4; \
  1662. mkl_dc_type tmp; \
  1663. MKL_INT jp; \
  1664. mkl_dc_type ajj; \
  1665. if (MKL_DC_IS_ZERO(a00)) { \
  1666. *info = 1; \
  1667. break; \
  1668. } \
  1669. if (MKL_DC_IS_ZERO(a11)) { \
  1670. *info = 2; \
  1671. break; \
  1672. } \
  1673. if (MKL_DC_IS_ZERO(a22)) { \
  1674. *info = 3; \
  1675. break; \
  1676. } \
  1677. if (MKL_DC_IS_ZERO(a33)) { \
  1678. *info = 4; \
  1679. break; \
  1680. } \
  1681. if (MKL_DC_IS_ZERO(a44)) { \
  1682. *info = 5; \
  1683. break; \
  1684. } \
  1685. *info = 0; \
  1686. MKL_DC_INV(a00, a00); \
  1687. MKL_DC_NEG(ajj, a00); \
  1688. MKL_DC_INV(a11, a11); \
  1689. MKL_DC_NEG(ajj, a11); \
  1690. tmp = a01; \
  1691. MKL_DC_MUL(a01, a01, a00); \
  1692. MKL_DC_MUL(a01, a01, ajj); \
  1693. MKL_DC_INV(a22, a22); \
  1694. MKL_DC_NEG(ajj, a22); \
  1695. tmp = a02; \
  1696. MKL_DC_MUL(a02, a02, a00); \
  1697. tmp = a12; \
  1698. MKL_DC_MUL_ADD(a02, tmp, a01, a02); \
  1699. MKL_DC_MUL(a12, a12, a11); \
  1700. MKL_DC_MUL(a02, a02, ajj); \
  1701. MKL_DC_MUL(a12, a12, ajj); \
  1702. MKL_DC_INV(a33, a33); \
  1703. MKL_DC_NEG(ajj, a33); \
  1704. tmp = a03; \
  1705. MKL_DC_MUL(a03, a03, a00); \
  1706. tmp = a13; \
  1707. MKL_DC_MUL_ADD(a03, tmp, a01, a03); \
  1708. MKL_DC_MUL(a13, a13, a11); \
  1709. tmp = a23; \
  1710. MKL_DC_MUL_ADD(a03, tmp, a02, a03); \
  1711. MKL_DC_MUL_ADD(a13, tmp, a12, a13); \
  1712. MKL_DC_MUL(a23, a23, a22); \
  1713. MKL_DC_MUL(a03, a03, ajj); \
  1714. MKL_DC_MUL(a13, a13, ajj); \
  1715. MKL_DC_MUL(a23, a23, ajj); \
  1716. MKL_DC_INV(a44, a44); \
  1717. MKL_DC_NEG(ajj, a44); \
  1718. tmp = a04; \
  1719. MKL_DC_MUL(a04, a04, a00); \
  1720. tmp = a14; \
  1721. MKL_DC_MUL_ADD(a04, tmp, a01, a04); \
  1722. MKL_DC_MUL(a14, a14, a11); \
  1723. tmp = a24; \
  1724. MKL_DC_MUL_ADD(a04, tmp, a02, a04); \
  1725. MKL_DC_MUL_ADD(a14, tmp, a12, a14); \
  1726. MKL_DC_MUL(a24, a24, a22); \
  1727. tmp = a34; \
  1728. MKL_DC_MUL_ADD(a04, tmp, a03, a04); \
  1729. MKL_DC_MUL_ADD(a14, tmp, a13, a14); \
  1730. MKL_DC_MUL_ADD(a24, tmp, a23, a24); \
  1731. MKL_DC_MUL(a34, a34, a33); \
  1732. MKL_DC_MUL(a04, a04, ajj); \
  1733. MKL_DC_MUL(a14, a14, ajj); \
  1734. MKL_DC_MUL(a24, a24, ajj); \
  1735. MKL_DC_MUL(a34, a34, ajj); \
  1736. w4 = a43; \
  1737. MKL_DC_SET_ZERO(a43); \
  1738. MKL_DC_NEG(tmp, w4); \
  1739. MKL_DC_MUL_ADD(a03, tmp, a04, a03); \
  1740. MKL_DC_MUL_ADD(a13, tmp, a14, a13); \
  1741. MKL_DC_MUL_ADD(a23, tmp, a24, a23); \
  1742. MKL_DC_MUL_ADD(a33, tmp, a34, a33); \
  1743. MKL_DC_MUL_ADD(a43, tmp, a44, a43); \
  1744. w3 = a32; \
  1745. MKL_DC_SET_ZERO(a32); \
  1746. w4 = a42; \
  1747. MKL_DC_SET_ZERO(a42); \
  1748. MKL_DC_NEG(tmp, w3); \
  1749. MKL_DC_MUL_ADD(a02, tmp, a03, a02); \
  1750. MKL_DC_MUL_ADD(a12, tmp, a13, a12); \
  1751. MKL_DC_MUL_ADD(a22, tmp, a23, a22); \
  1752. MKL_DC_MUL_ADD(a32, tmp, a33, a32); \
  1753. MKL_DC_MUL_ADD(a42, tmp, a43, a42); \
  1754. MKL_DC_NEG(tmp, w4); \
  1755. MKL_DC_MUL_ADD(a02, tmp, a04, a02); \
  1756. MKL_DC_MUL_ADD(a12, tmp, a14, a12); \
  1757. MKL_DC_MUL_ADD(a22, tmp, a24, a22); \
  1758. MKL_DC_MUL_ADD(a32, tmp, a34, a32); \
  1759. MKL_DC_MUL_ADD(a42, tmp, a44, a42); \
  1760. w2 = a21; \
  1761. MKL_DC_SET_ZERO(a21); \
  1762. w3 = a31; \
  1763. MKL_DC_SET_ZERO(a31); \
  1764. w4 = a41; \
  1765. MKL_DC_SET_ZERO(a41); \
  1766. MKL_DC_NEG(tmp, w2); \
  1767. MKL_DC_MUL_ADD(a01, tmp, a02, a01); \
  1768. MKL_DC_MUL_ADD(a11, tmp, a12, a11); \
  1769. MKL_DC_MUL_ADD(a21, tmp, a22, a21); \
  1770. MKL_DC_MUL_ADD(a31, tmp, a32, a31); \
  1771. MKL_DC_MUL_ADD(a41, tmp, a42, a41); \
  1772. MKL_DC_NEG(tmp, w3); \
  1773. MKL_DC_MUL_ADD(a01, tmp, a03, a01); \
  1774. MKL_DC_MUL_ADD(a11, tmp, a13, a11); \
  1775. MKL_DC_MUL_ADD(a21, tmp, a23, a21); \
  1776. MKL_DC_MUL_ADD(a31, tmp, a33, a31); \
  1777. MKL_DC_MUL_ADD(a41, tmp, a43, a41); \
  1778. MKL_DC_NEG(tmp, w4); \
  1779. MKL_DC_MUL_ADD(a01, tmp, a04, a01); \
  1780. MKL_DC_MUL_ADD(a11, tmp, a14, a11); \
  1781. MKL_DC_MUL_ADD(a21, tmp, a24, a21); \
  1782. MKL_DC_MUL_ADD(a31, tmp, a34, a31); \
  1783. MKL_DC_MUL_ADD(a41, tmp, a44, a41); \
  1784. w1 = a10; \
  1785. MKL_DC_SET_ZERO(a10); \
  1786. w2 = a20; \
  1787. MKL_DC_SET_ZERO(a20); \
  1788. w3 = a30; \
  1789. MKL_DC_SET_ZERO(a30); \
  1790. w4 = a40; \
  1791. MKL_DC_SET_ZERO(a40); \
  1792. MKL_DC_NEG(tmp, w1); \
  1793. MKL_DC_MUL_ADD(a00, tmp, a01, a00); \
  1794. MKL_DC_MUL_ADD(a10, tmp, a11, a10); \
  1795. MKL_DC_MUL_ADD(a20, tmp, a21, a20); \
  1796. MKL_DC_MUL_ADD(a30, tmp, a31, a30); \
  1797. MKL_DC_MUL_ADD(a40, tmp, a41, a40); \
  1798. MKL_DC_NEG(tmp, w2); \
  1799. MKL_DC_MUL_ADD(a00, tmp, a02, a00); \
  1800. MKL_DC_MUL_ADD(a10, tmp, a12, a10); \
  1801. MKL_DC_MUL_ADD(a20, tmp, a22, a20); \
  1802. MKL_DC_MUL_ADD(a30, tmp, a32, a30); \
  1803. MKL_DC_MUL_ADD(a40, tmp, a42, a40); \
  1804. MKL_DC_NEG(tmp, w3); \
  1805. MKL_DC_MUL_ADD(a00, tmp, a03, a00); \
  1806. MKL_DC_MUL_ADD(a10, tmp, a13, a10); \
  1807. MKL_DC_MUL_ADD(a20, tmp, a23, a20); \
  1808. MKL_DC_MUL_ADD(a30, tmp, a33, a30); \
  1809. MKL_DC_MUL_ADD(a40, tmp, a43, a40); \
  1810. MKL_DC_NEG(tmp, w4); \
  1811. MKL_DC_MUL_ADD(a00, tmp, a04, a00); \
  1812. MKL_DC_MUL_ADD(a10, tmp, a14, a10); \
  1813. MKL_DC_MUL_ADD(a20, tmp, a24, a20); \
  1814. MKL_DC_MUL_ADD(a30, tmp, a34, a30); \
  1815. MKL_DC_MUL_ADD(a40, tmp, a44, a40); \
  1816. dc_access(a, lda, 0, 0) = a00; \
  1817. dc_access(a, lda, 0, 1) = a01; \
  1818. dc_access(a, lda, 0, 2) = a02; \
  1819. dc_access(a, lda, 0, 3) = a03; \
  1820. dc_access(a, lda, 0, 4) = a04; \
  1821. dc_access(a, lda, 1, 0) = a10; \
  1822. dc_access(a, lda, 1, 1) = a11; \
  1823. dc_access(a, lda, 1, 2) = a12; \
  1824. dc_access(a, lda, 1, 3) = a13; \
  1825. dc_access(a, lda, 1, 4) = a14; \
  1826. dc_access(a, lda, 2, 0) = a20; \
  1827. dc_access(a, lda, 2, 1) = a21; \
  1828. dc_access(a, lda, 2, 2) = a22; \
  1829. dc_access(a, lda, 2, 3) = a23; \
  1830. dc_access(a, lda, 2, 4) = a24; \
  1831. dc_access(a, lda, 3, 0) = a30; \
  1832. dc_access(a, lda, 3, 1) = a31; \
  1833. dc_access(a, lda, 3, 2) = a32; \
  1834. dc_access(a, lda, 3, 3) = a33; \
  1835. dc_access(a, lda, 3, 4) = a34; \
  1836. dc_access(a, lda, 4, 0) = a40; \
  1837. dc_access(a, lda, 4, 1) = a41; \
  1838. dc_access(a, lda, 4, 2) = a42; \
  1839. dc_access(a, lda, 4, 3) = a43; \
  1840. dc_access(a, lda, 4, 4) = a44; \
  1841. jp = ipiv[3] - 1; \
  1842. if (jp != 3) { \
  1843. MKL_DC_SWAP(dc_access(a, lda, 0, 3), dc_access(a, lda, 0, jp)); \
  1844. MKL_DC_SWAP(dc_access(a, lda, 1, 3), dc_access(a, lda, 1, jp)); \
  1845. MKL_DC_SWAP(dc_access(a, lda, 2, 3), dc_access(a, lda, 2, jp)); \
  1846. MKL_DC_SWAP(dc_access(a, lda, 3, 3), dc_access(a, lda, 3, jp)); \
  1847. MKL_DC_SWAP(dc_access(a, lda, 4, 3), dc_access(a, lda, 4, jp)); \
  1848. } \
  1849. jp = ipiv[2] - 1; \
  1850. if (jp != 2) { \
  1851. MKL_DC_SWAP(dc_access(a, lda, 0, 2), dc_access(a, lda, 0, jp)); \
  1852. MKL_DC_SWAP(dc_access(a, lda, 1, 2), dc_access(a, lda, 1, jp)); \
  1853. MKL_DC_SWAP(dc_access(a, lda, 2, 2), dc_access(a, lda, 2, jp)); \
  1854. MKL_DC_SWAP(dc_access(a, lda, 3, 2), dc_access(a, lda, 3, jp)); \
  1855. MKL_DC_SWAP(dc_access(a, lda, 4, 2), dc_access(a, lda, 4, jp)); \
  1856. } \
  1857. jp = ipiv[1] - 1; \
  1858. if (jp != 1) { \
  1859. MKL_DC_SWAP(dc_access(a, lda, 0, 1), dc_access(a, lda, 0, jp)); \
  1860. MKL_DC_SWAP(dc_access(a, lda, 1, 1), dc_access(a, lda, 1, jp)); \
  1861. MKL_DC_SWAP(dc_access(a, lda, 2, 1), dc_access(a, lda, 2, jp)); \
  1862. MKL_DC_SWAP(dc_access(a, lda, 3, 1), dc_access(a, lda, 3, jp)); \
  1863. MKL_DC_SWAP(dc_access(a, lda, 4, 1), dc_access(a, lda, 4, jp)); \
  1864. } \
  1865. jp = ipiv[0] - 1; \
  1866. if (jp != 0) { \
  1867. MKL_DC_SWAP(dc_access(a, lda, 0, 0), dc_access(a, lda, 0, jp)); \
  1868. MKL_DC_SWAP(dc_access(a, lda, 1, 0), dc_access(a, lda, 1, jp)); \
  1869. MKL_DC_SWAP(dc_access(a, lda, 2, 0), dc_access(a, lda, 2, jp)); \
  1870. MKL_DC_SWAP(dc_access(a, lda, 3, 0), dc_access(a, lda, 3, jp)); \
  1871. MKL_DC_SWAP(dc_access(a, lda, 4, 0), dc_access(a, lda, 4, jp)); \
  1872. } \
  1873. } while (0)
  1874. #define mkl_dc_getri_6x6(n, a, lda, ipiv, info, dc_access) \
  1875. do { \
  1876. mkl_dc_type a00 = dc_access(a, lda, 0, 0); \
  1877. mkl_dc_type a01 = dc_access(a, lda, 0, 1); \
  1878. mkl_dc_type a02 = dc_access(a, lda, 0, 2); \
  1879. mkl_dc_type a03 = dc_access(a, lda, 0, 3); \
  1880. mkl_dc_type a04 = dc_access(a, lda, 0, 4); \
  1881. mkl_dc_type a05 = dc_access(a, lda, 0, 5); \
  1882. mkl_dc_type a10 = dc_access(a, lda, 1, 0); \
  1883. mkl_dc_type a11 = dc_access(a, lda, 1, 1); \
  1884. mkl_dc_type a12 = dc_access(a, lda, 1, 2); \
  1885. mkl_dc_type a13 = dc_access(a, lda, 1, 3); \
  1886. mkl_dc_type a14 = dc_access(a, lda, 1, 4); \
  1887. mkl_dc_type a15 = dc_access(a, lda, 1, 5); \
  1888. mkl_dc_type a20 = dc_access(a, lda, 2, 0); \
  1889. mkl_dc_type a21 = dc_access(a, lda, 2, 1); \
  1890. mkl_dc_type a22 = dc_access(a, lda, 2, 2); \
  1891. mkl_dc_type a23 = dc_access(a, lda, 2, 3); \
  1892. mkl_dc_type a24 = dc_access(a, lda, 2, 4); \
  1893. mkl_dc_type a25 = dc_access(a, lda, 2, 5); \
  1894. mkl_dc_type a30 = dc_access(a, lda, 3, 0); \
  1895. mkl_dc_type a31 = dc_access(a, lda, 3, 1); \
  1896. mkl_dc_type a32 = dc_access(a, lda, 3, 2); \
  1897. mkl_dc_type a33 = dc_access(a, lda, 3, 3); \
  1898. mkl_dc_type a34 = dc_access(a, lda, 3, 4); \
  1899. mkl_dc_type a35 = dc_access(a, lda, 3, 5); \
  1900. mkl_dc_type a40 = dc_access(a, lda, 4, 0); \
  1901. mkl_dc_type a41 = dc_access(a, lda, 4, 1); \
  1902. mkl_dc_type a42 = dc_access(a, lda, 4, 2); \
  1903. mkl_dc_type a43 = dc_access(a, lda, 4, 3); \
  1904. mkl_dc_type a44 = dc_access(a, lda, 4, 4); \
  1905. mkl_dc_type a45 = dc_access(a, lda, 4, 5); \
  1906. mkl_dc_type a50 = dc_access(a, lda, 5, 0); \
  1907. mkl_dc_type a51 = dc_access(a, lda, 5, 1); \
  1908. mkl_dc_type a52 = dc_access(a, lda, 5, 2); \
  1909. mkl_dc_type a53 = dc_access(a, lda, 5, 3); \
  1910. mkl_dc_type a54 = dc_access(a, lda, 5, 4); \
  1911. mkl_dc_type a55 = dc_access(a, lda, 5, 5); \
  1912. mkl_dc_type w1; \
  1913. mkl_dc_type w2; \
  1914. mkl_dc_type w3; \
  1915. mkl_dc_type w4; \
  1916. mkl_dc_type w5; \
  1917. mkl_dc_type tmp; \
  1918. MKL_INT jp; \
  1919. mkl_dc_type ajj; \
  1920. if (MKL_DC_IS_ZERO(a00)) { \
  1921. *info = 1; \
  1922. break; \
  1923. } \
  1924. if (MKL_DC_IS_ZERO(a11)) { \
  1925. *info = 2; \
  1926. break; \
  1927. } \
  1928. if (MKL_DC_IS_ZERO(a22)) { \
  1929. *info = 3; \
  1930. break; \
  1931. } \
  1932. if (MKL_DC_IS_ZERO(a33)) { \
  1933. *info = 4; \
  1934. break; \
  1935. } \
  1936. if (MKL_DC_IS_ZERO(a44)) { \
  1937. *info = 5; \
  1938. break; \
  1939. } \
  1940. if (MKL_DC_IS_ZERO(a55)) { \
  1941. *info = 6; \
  1942. break; \
  1943. } \
  1944. *info = 0; \
  1945. MKL_DC_INV(a00, a00); \
  1946. MKL_DC_NEG(ajj, a00); \
  1947. MKL_DC_INV(a11, a11); \
  1948. MKL_DC_NEG(ajj, a11); \
  1949. tmp = a01; \
  1950. MKL_DC_MUL(a01, a01, a00); \
  1951. MKL_DC_MUL(a01, a01, ajj); \
  1952. MKL_DC_INV(a22, a22); \
  1953. MKL_DC_NEG(ajj, a22); \
  1954. tmp = a02; \
  1955. MKL_DC_MUL(a02, a02, a00); \
  1956. tmp = a12; \
  1957. MKL_DC_MUL_ADD(a02, tmp, a01, a02); \
  1958. MKL_DC_MUL(a12, a12, a11); \
  1959. MKL_DC_MUL(a02, a02, ajj); \
  1960. MKL_DC_MUL(a12, a12, ajj); \
  1961. MKL_DC_INV(a33, a33); \
  1962. MKL_DC_NEG(ajj, a33); \
  1963. tmp = a03; \
  1964. MKL_DC_MUL(a03, a03, a00); \
  1965. tmp = a13; \
  1966. MKL_DC_MUL_ADD(a03, tmp, a01, a03); \
  1967. MKL_DC_MUL(a13, a13, a11); \
  1968. tmp = a23; \
  1969. MKL_DC_MUL_ADD(a03, tmp, a02, a03); \
  1970. MKL_DC_MUL_ADD(a13, tmp, a12, a13); \
  1971. MKL_DC_MUL(a23, a23, a22); \
  1972. MKL_DC_MUL(a03, a03, ajj); \
  1973. MKL_DC_MUL(a13, a13, ajj); \
  1974. MKL_DC_MUL(a23, a23, ajj); \
  1975. MKL_DC_INV(a44, a44); \
  1976. MKL_DC_NEG(ajj, a44); \
  1977. tmp = a04; \
  1978. MKL_DC_MUL(a04, a04, a00); \
  1979. tmp = a14; \
  1980. MKL_DC_MUL_ADD(a04, tmp, a01, a04); \
  1981. MKL_DC_MUL(a14, a14, a11); \
  1982. tmp = a24; \
  1983. MKL_DC_MUL_ADD(a04, tmp, a02, a04); \
  1984. MKL_DC_MUL_ADD(a14, tmp, a12, a14); \
  1985. MKL_DC_MUL(a24, a24, a22); \
  1986. tmp = a34; \
  1987. MKL_DC_MUL_ADD(a04, tmp, a03, a04); \
  1988. MKL_DC_MUL_ADD(a14, tmp, a13, a14); \
  1989. MKL_DC_MUL_ADD(a24, tmp, a23, a24); \
  1990. MKL_DC_MUL(a34, a34, a33); \
  1991. MKL_DC_MUL(a04, a04, ajj); \
  1992. MKL_DC_MUL(a14, a14, ajj); \
  1993. MKL_DC_MUL(a24, a24, ajj); \
  1994. MKL_DC_MUL(a34, a34, ajj); \
  1995. MKL_DC_INV(a55, a55); \
  1996. MKL_DC_NEG(ajj, a55); \
  1997. tmp = a05; \
  1998. MKL_DC_MUL(a05, a05, a00); \
  1999. tmp = a15; \
  2000. MKL_DC_MUL_ADD(a05, tmp, a01, a05); \
  2001. MKL_DC_MUL(a15, a15, a11); \
  2002. tmp = a25; \
  2003. MKL_DC_MUL_ADD(a05, tmp, a02, a05); \
  2004. MKL_DC_MUL_ADD(a15, tmp, a12, a15); \
  2005. MKL_DC_MUL(a25, a25, a22); \
  2006. tmp = a35; \
  2007. MKL_DC_MUL_ADD(a05, tmp, a03, a05); \
  2008. MKL_DC_MUL_ADD(a15, tmp, a13, a15); \
  2009. MKL_DC_MUL_ADD(a25, tmp, a23, a25); \
  2010. MKL_DC_MUL(a35, a35, a33); \
  2011. tmp = a45; \
  2012. MKL_DC_MUL_ADD(a05, tmp, a04, a05); \
  2013. MKL_DC_MUL_ADD(a15, tmp, a14, a15); \
  2014. MKL_DC_MUL_ADD(a25, tmp, a24, a25); \
  2015. MKL_DC_MUL_ADD(a35, tmp, a34, a35); \
  2016. MKL_DC_MUL(a45, a45, a44); \
  2017. MKL_DC_MUL(a05, a05, ajj); \
  2018. MKL_DC_MUL(a15, a15, ajj); \
  2019. MKL_DC_MUL(a25, a25, ajj); \
  2020. MKL_DC_MUL(a35, a35, ajj); \
  2021. MKL_DC_MUL(a45, a45, ajj); \
  2022. w5 = a54; \
  2023. MKL_DC_SET_ZERO(a54); \
  2024. MKL_DC_NEG(tmp, w5); \
  2025. MKL_DC_MUL_ADD(a04, tmp, a05, a04); \
  2026. MKL_DC_MUL_ADD(a14, tmp, a15, a14); \
  2027. MKL_DC_MUL_ADD(a24, tmp, a25, a24); \
  2028. MKL_DC_MUL_ADD(a34, tmp, a35, a34); \
  2029. MKL_DC_MUL_ADD(a44, tmp, a45, a44); \
  2030. MKL_DC_MUL_ADD(a54, tmp, a55, a54); \
  2031. w4 = a43; \
  2032. MKL_DC_SET_ZERO(a43); \
  2033. w5 = a53; \
  2034. MKL_DC_SET_ZERO(a53); \
  2035. MKL_DC_NEG(tmp, w4); \
  2036. MKL_DC_MUL_ADD(a03, tmp, a04, a03); \
  2037. MKL_DC_MUL_ADD(a13, tmp, a14, a13); \
  2038. MKL_DC_MUL_ADD(a23, tmp, a24, a23); \
  2039. MKL_DC_MUL_ADD(a33, tmp, a34, a33); \
  2040. MKL_DC_MUL_ADD(a43, tmp, a44, a43); \
  2041. MKL_DC_MUL_ADD(a53, tmp, a54, a53); \
  2042. MKL_DC_NEG(tmp, w5); \
  2043. MKL_DC_MUL_ADD(a03, tmp, a05, a03); \
  2044. MKL_DC_MUL_ADD(a13, tmp, a15, a13); \
  2045. MKL_DC_MUL_ADD(a23, tmp, a25, a23); \
  2046. MKL_DC_MUL_ADD(a33, tmp, a35, a33); \
  2047. MKL_DC_MUL_ADD(a43, tmp, a45, a43); \
  2048. MKL_DC_MUL_ADD(a53, tmp, a55, a53); \
  2049. w3 = a32; \
  2050. MKL_DC_SET_ZERO(a32); \
  2051. w4 = a42; \
  2052. MKL_DC_SET_ZERO(a42); \
  2053. w5 = a52; \
  2054. MKL_DC_SET_ZERO(a52); \
  2055. MKL_DC_NEG(tmp, w3); \
  2056. MKL_DC_MUL_ADD(a02, tmp, a03, a02); \
  2057. MKL_DC_MUL_ADD(a12, tmp, a13, a12); \
  2058. MKL_DC_MUL_ADD(a22, tmp, a23, a22); \
  2059. MKL_DC_MUL_ADD(a32, tmp, a33, a32); \
  2060. MKL_DC_MUL_ADD(a42, tmp, a43, a42); \
  2061. MKL_DC_MUL_ADD(a52, tmp, a53, a52); \
  2062. MKL_DC_NEG(tmp, w4); \
  2063. MKL_DC_MUL_ADD(a02, tmp, a04, a02); \
  2064. MKL_DC_MUL_ADD(a12, tmp, a14, a12); \
  2065. MKL_DC_MUL_ADD(a22, tmp, a24, a22); \
  2066. MKL_DC_MUL_ADD(a32, tmp, a34, a32); \
  2067. MKL_DC_MUL_ADD(a42, tmp, a44, a42); \
  2068. MKL_DC_MUL_ADD(a52, tmp, a54, a52); \
  2069. MKL_DC_NEG(tmp, w5); \
  2070. MKL_DC_MUL_ADD(a02, tmp, a05, a02); \
  2071. MKL_DC_MUL_ADD(a12, tmp, a15, a12); \
  2072. MKL_DC_MUL_ADD(a22, tmp, a25, a22); \
  2073. MKL_DC_MUL_ADD(a32, tmp, a35, a32); \
  2074. MKL_DC_MUL_ADD(a42, tmp, a45, a42); \
  2075. MKL_DC_MUL_ADD(a52, tmp, a55, a52); \
  2076. w2 = a21; \
  2077. MKL_DC_SET_ZERO(a21); \
  2078. w3 = a31; \
  2079. MKL_DC_SET_ZERO(a31); \
  2080. w4 = a41; \
  2081. MKL_DC_SET_ZERO(a41); \
  2082. w5 = a51; \
  2083. MKL_DC_SET_ZERO(a51); \
  2084. MKL_DC_NEG(tmp, w2); \
  2085. MKL_DC_MUL_ADD(a01, tmp, a02, a01); \
  2086. MKL_DC_MUL_ADD(a11, tmp, a12, a11); \
  2087. MKL_DC_MUL_ADD(a21, tmp, a22, a21); \
  2088. MKL_DC_MUL_ADD(a31, tmp, a32, a31); \
  2089. MKL_DC_MUL_ADD(a41, tmp, a42, a41); \
  2090. MKL_DC_MUL_ADD(a51, tmp, a52, a51); \
  2091. MKL_DC_NEG(tmp, w3); \
  2092. MKL_DC_MUL_ADD(a01, tmp, a03, a01); \
  2093. MKL_DC_MUL_ADD(a11, tmp, a13, a11); \
  2094. MKL_DC_MUL_ADD(a21, tmp, a23, a21); \
  2095. MKL_DC_MUL_ADD(a31, tmp, a33, a31); \
  2096. MKL_DC_MUL_ADD(a41, tmp, a43, a41); \
  2097. MKL_DC_MUL_ADD(a51, tmp, a53, a51); \
  2098. MKL_DC_NEG(tmp, w4); \
  2099. MKL_DC_MUL_ADD(a01, tmp, a04, a01); \
  2100. MKL_DC_MUL_ADD(a11, tmp, a14, a11); \
  2101. MKL_DC_MUL_ADD(a21, tmp, a24, a21); \
  2102. MKL_DC_MUL_ADD(a31, tmp, a34, a31); \
  2103. MKL_DC_MUL_ADD(a41, tmp, a44, a41); \
  2104. MKL_DC_MUL_ADD(a51, tmp, a54, a51); \
  2105. MKL_DC_NEG(tmp, w5); \
  2106. MKL_DC_MUL_ADD(a01, tmp, a05, a01); \
  2107. MKL_DC_MUL_ADD(a11, tmp, a15, a11); \
  2108. MKL_DC_MUL_ADD(a21, tmp, a25, a21); \
  2109. MKL_DC_MUL_ADD(a31, tmp, a35, a31); \
  2110. MKL_DC_MUL_ADD(a41, tmp, a45, a41); \
  2111. MKL_DC_MUL_ADD(a51, tmp, a55, a51); \
  2112. w1 = a10; \
  2113. MKL_DC_SET_ZERO(a10); \
  2114. w2 = a20; \
  2115. MKL_DC_SET_ZERO(a20); \
  2116. w3 = a30; \
  2117. MKL_DC_SET_ZERO(a30); \
  2118. w4 = a40; \
  2119. MKL_DC_SET_ZERO(a40); \
  2120. w5 = a50; \
  2121. MKL_DC_SET_ZERO(a50); \
  2122. MKL_DC_NEG(tmp, w1); \
  2123. MKL_DC_MUL_ADD(a00, tmp, a01, a00); \
  2124. MKL_DC_MUL_ADD(a10, tmp, a11, a10); \
  2125. MKL_DC_MUL_ADD(a20, tmp, a21, a20); \
  2126. MKL_DC_MUL_ADD(a30, tmp, a31, a30); \
  2127. MKL_DC_MUL_ADD(a40, tmp, a41, a40); \
  2128. MKL_DC_MUL_ADD(a50, tmp, a51, a50); \
  2129. MKL_DC_NEG(tmp, w2); \
  2130. MKL_DC_MUL_ADD(a00, tmp, a02, a00); \
  2131. MKL_DC_MUL_ADD(a10, tmp, a12, a10); \
  2132. MKL_DC_MUL_ADD(a20, tmp, a22, a20); \
  2133. MKL_DC_MUL_ADD(a30, tmp, a32, a30); \
  2134. MKL_DC_MUL_ADD(a40, tmp, a42, a40); \
  2135. MKL_DC_MUL_ADD(a50, tmp, a52, a50); \
  2136. MKL_DC_NEG(tmp, w3); \
  2137. MKL_DC_MUL_ADD(a00, tmp, a03, a00); \
  2138. MKL_DC_MUL_ADD(a10, tmp, a13, a10); \
  2139. MKL_DC_MUL_ADD(a20, tmp, a23, a20); \
  2140. MKL_DC_MUL_ADD(a30, tmp, a33, a30); \
  2141. MKL_DC_MUL_ADD(a40, tmp, a43, a40); \
  2142. MKL_DC_MUL_ADD(a50, tmp, a53, a50); \
  2143. MKL_DC_NEG(tmp, w4); \
  2144. MKL_DC_MUL_ADD(a00, tmp, a04, a00); \
  2145. MKL_DC_MUL_ADD(a10, tmp, a14, a10); \
  2146. MKL_DC_MUL_ADD(a20, tmp, a24, a20); \
  2147. MKL_DC_MUL_ADD(a30, tmp, a34, a30); \
  2148. MKL_DC_MUL_ADD(a40, tmp, a44, a40); \
  2149. MKL_DC_MUL_ADD(a50, tmp, a54, a50); \
  2150. MKL_DC_NEG(tmp, w5); \
  2151. MKL_DC_MUL_ADD(a00, tmp, a05, a00); \
  2152. MKL_DC_MUL_ADD(a10, tmp, a15, a10); \
  2153. MKL_DC_MUL_ADD(a20, tmp, a25, a20); \
  2154. MKL_DC_MUL_ADD(a30, tmp, a35, a30); \
  2155. MKL_DC_MUL_ADD(a40, tmp, a45, a40); \
  2156. MKL_DC_MUL_ADD(a50, tmp, a55, a50); \
  2157. dc_access(a, lda, 0, 0) = a00; \
  2158. dc_access(a, lda, 0, 1) = a01; \
  2159. dc_access(a, lda, 0, 2) = a02; \
  2160. dc_access(a, lda, 0, 3) = a03; \
  2161. dc_access(a, lda, 0, 4) = a04; \
  2162. dc_access(a, lda, 0, 5) = a05; \
  2163. dc_access(a, lda, 1, 0) = a10; \
  2164. dc_access(a, lda, 1, 1) = a11; \
  2165. dc_access(a, lda, 1, 2) = a12; \
  2166. dc_access(a, lda, 1, 3) = a13; \
  2167. dc_access(a, lda, 1, 4) = a14; \
  2168. dc_access(a, lda, 1, 5) = a15; \
  2169. dc_access(a, lda, 2, 0) = a20; \
  2170. dc_access(a, lda, 2, 1) = a21; \
  2171. dc_access(a, lda, 2, 2) = a22; \
  2172. dc_access(a, lda, 2, 3) = a23; \
  2173. dc_access(a, lda, 2, 4) = a24; \
  2174. dc_access(a, lda, 2, 5) = a25; \
  2175. dc_access(a, lda, 3, 0) = a30; \
  2176. dc_access(a, lda, 3, 1) = a31; \
  2177. dc_access(a, lda, 3, 2) = a32; \
  2178. dc_access(a, lda, 3, 3) = a33; \
  2179. dc_access(a, lda, 3, 4) = a34; \
  2180. dc_access(a, lda, 3, 5) = a35; \
  2181. dc_access(a, lda, 4, 0) = a40; \
  2182. dc_access(a, lda, 4, 1) = a41; \
  2183. dc_access(a, lda, 4, 2) = a42; \
  2184. dc_access(a, lda, 4, 3) = a43; \
  2185. dc_access(a, lda, 4, 4) = a44; \
  2186. dc_access(a, lda, 4, 5) = a45; \
  2187. dc_access(a, lda, 5, 0) = a50; \
  2188. dc_access(a, lda, 5, 1) = a51; \
  2189. dc_access(a, lda, 5, 2) = a52; \
  2190. dc_access(a, lda, 5, 3) = a53; \
  2191. dc_access(a, lda, 5, 4) = a54; \
  2192. dc_access(a, lda, 5, 5) = a55; \
  2193. jp = ipiv[4] - 1; \
  2194. if (jp != 4) { \
  2195. MKL_DC_SWAP(dc_access(a, lda, 0, 4), dc_access(a, lda, 0, jp)); \
  2196. MKL_DC_SWAP(dc_access(a, lda, 1, 4), dc_access(a, lda, 1, jp)); \
  2197. MKL_DC_SWAP(dc_access(a, lda, 2, 4), dc_access(a, lda, 2, jp)); \
  2198. MKL_DC_SWAP(dc_access(a, lda, 3, 4), dc_access(a, lda, 3, jp)); \
  2199. MKL_DC_SWAP(dc_access(a, lda, 4, 4), dc_access(a, lda, 4, jp)); \
  2200. MKL_DC_SWAP(dc_access(a, lda, 5, 4), dc_access(a, lda, 5, jp)); \
  2201. } \
  2202. jp = ipiv[3] - 1; \
  2203. if (jp != 3) { \
  2204. MKL_DC_SWAP(dc_access(a, lda, 0, 3), dc_access(a, lda, 0, jp)); \
  2205. MKL_DC_SWAP(dc_access(a, lda, 1, 3), dc_access(a, lda, 1, jp)); \
  2206. MKL_DC_SWAP(dc_access(a, lda, 2, 3), dc_access(a, lda, 2, jp)); \
  2207. MKL_DC_SWAP(dc_access(a, lda, 3, 3), dc_access(a, lda, 3, jp)); \
  2208. MKL_DC_SWAP(dc_access(a, lda, 4, 3), dc_access(a, lda, 4, jp)); \
  2209. MKL_DC_SWAP(dc_access(a, lda, 5, 3), dc_access(a, lda, 5, jp)); \
  2210. } \
  2211. jp = ipiv[2] - 1; \
  2212. if (jp != 2) { \
  2213. MKL_DC_SWAP(dc_access(a, lda, 0, 2), dc_access(a, lda, 0, jp)); \
  2214. MKL_DC_SWAP(dc_access(a, lda, 1, 2), dc_access(a, lda, 1, jp)); \
  2215. MKL_DC_SWAP(dc_access(a, lda, 2, 2), dc_access(a, lda, 2, jp)); \
  2216. MKL_DC_SWAP(dc_access(a, lda, 3, 2), dc_access(a, lda, 3, jp)); \
  2217. MKL_DC_SWAP(dc_access(a, lda, 4, 2), dc_access(a, lda, 4, jp)); \
  2218. MKL_DC_SWAP(dc_access(a, lda, 5, 2), dc_access(a, lda, 5, jp)); \
  2219. } \
  2220. jp = ipiv[1] - 1; \
  2221. if (jp != 1) { \
  2222. MKL_DC_SWAP(dc_access(a, lda, 0, 1), dc_access(a, lda, 0, jp)); \
  2223. MKL_DC_SWAP(dc_access(a, lda, 1, 1), dc_access(a, lda, 1, jp)); \
  2224. MKL_DC_SWAP(dc_access(a, lda, 2, 1), dc_access(a, lda, 2, jp)); \
  2225. MKL_DC_SWAP(dc_access(a, lda, 3, 1), dc_access(a, lda, 3, jp)); \
  2226. MKL_DC_SWAP(dc_access(a, lda, 4, 1), dc_access(a, lda, 4, jp)); \
  2227. MKL_DC_SWAP(dc_access(a, lda, 5, 1), dc_access(a, lda, 5, jp)); \
  2228. } \
  2229. jp = ipiv[0] - 1; \
  2230. if (jp != 0) { \
  2231. MKL_DC_SWAP(dc_access(a, lda, 0, 0), dc_access(a, lda, 0, jp)); \
  2232. MKL_DC_SWAP(dc_access(a, lda, 1, 0), dc_access(a, lda, 1, jp)); \
  2233. MKL_DC_SWAP(dc_access(a, lda, 2, 0), dc_access(a, lda, 2, jp)); \
  2234. MKL_DC_SWAP(dc_access(a, lda, 3, 0), dc_access(a, lda, 3, jp)); \
  2235. MKL_DC_SWAP(dc_access(a, lda, 4, 0), dc_access(a, lda, 4, jp)); \
  2236. MKL_DC_SWAP(dc_access(a, lda, 5, 0), dc_access(a, lda, 5, jp)); \
  2237. } \
  2238. } while (0)
  2239. #define mkl_dc_getri_inverse_upper(n, a, lda, dc_access) \
  2240. do { \
  2241. MKL_INT i, j, ii, jj; \
  2242. mkl_dc_type ajj; \
  2243. for (j = 0; j < n; ++j) { \
  2244. MKL_DC_INV(dc_access(a, lda, j, j), dc_access(a, lda, j, j)); \
  2245. MKL_DC_NEG(ajj, dc_access(a, lda, j, j)); \
  2246. for (jj = 0; jj < j; ++jj) { \
  2247. mkl_dc_type tmp = dc_access(a, lda, jj, j); \
  2248. for (ii = 0; ii < jj; ++ii) { \
  2249. MKL_DC_MUL_ADD(dc_access(a, lda, ii, j), tmp, dc_access(a, lda, ii, jj), dc_access(a, lda, ii, j)); \
  2250. } \
  2251. MKL_DC_MUL(dc_access(a, lda, jj, j), dc_access(a, lda, jj, j), dc_access(a, lda, jj, jj)); \
  2252. } \
  2253. for (i = 0; i < j; ++i) { \
  2254. MKL_DC_MUL(dc_access(a, lda, i, j), dc_access(a, lda, i, j), ajj); \
  2255. } \
  2256. } \
  2257. } while (0)
  2258. #define mkl_dc_getri_solve_lower(n, a, lda, work, dc_access) \
  2259. do { \
  2260. MKL_INT i, j, ii, jj; \
  2261. mkl_dc_type tmp; \
  2262. for (j = n - 1; j >= 0; j--) { \
  2263. for (i = j + 1; i < n; ++i) { \
  2264. work[i] = dc_access(a, lda, i, j); \
  2265. MKL_DC_SET_ZERO(dc_access(a, lda, i, j)); \
  2266. } \
  2267. for (jj = j + 1; jj < n; ++jj) { \
  2268. MKL_DC_NEG(tmp, work[jj]); \
  2269. for (ii = 0; ii < n; ++ii) { \
  2270. MKL_DC_MUL_ADD(dc_access(a, lda, ii, j), tmp, dc_access(a, lda, ii, jj), dc_access(a, lda, ii, j)); \
  2271. } \
  2272. } \
  2273. } \
  2274. } while (0)
  2275. #define mkl_dc_getri_n(n, a, lda, ipiv, work, lwork, info, dc_access) \
  2276. do { \
  2277. MKL_INT i, j, jp; \
  2278. *info = 0; \
  2279. for (i = 0; i < n; ++i) { \
  2280. if (MKL_DC_IS_ZERO(dc_access(a, lda, i, i))) { \
  2281. *info = i + 1; \
  2282. break; \
  2283. } \
  2284. } \
  2285. if (*info == 0) { \
  2286. mkl_dc_getri_inverse_upper(n, a, lda, dc_access); \
  2287. mkl_dc_getri_solve_lower(n, a, lda, work, dc_access); \
  2288. for (j = n - 2; j >= 0; j--) { \
  2289. jp = ipiv[j] - 1; \
  2290. if (jp != j) { \
  2291. for (i = 0; i < n; ++i) { \
  2292. MKL_DC_SWAP(dc_access(a, lda, i, j), dc_access(a, lda, i, jp)); \
  2293. } \
  2294. } \
  2295. } \
  2296. } \
  2297. } while (0)
  2298. #define mkl_dc_getri_generic(n, a, lda, ipiv, work, lwork, info, dc_access) \
  2299. do { \
  2300. if (lwork != -1) { \
  2301. switch (n) { \
  2302. case 0: \
  2303. break; \
  2304. case 1: \
  2305. mkl_dc_getri_1x1(1, a, lda, ipiv, info, dc_access); \
  2306. break; \
  2307. case 2: \
  2308. mkl_dc_getri_2x2(2, a, lda, ipiv, info, dc_access); \
  2309. break; \
  2310. case 3: \
  2311. mkl_dc_getri_3x3(3, a, lda, ipiv, info, dc_access); \
  2312. break; \
  2313. case 4: \
  2314. mkl_dc_getri_4x4(4, a, lda, ipiv, info, dc_access); \
  2315. break; \
  2316. case 5: \
  2317. mkl_dc_getri_5x5(5, a, lda, ipiv, info, dc_access); \
  2318. break; \
  2319. case 6: \
  2320. mkl_dc_getri_6x6(6, a, lda, ipiv, info, dc_access); \
  2321. break; \
  2322. default: \
  2323. mkl_dc_getri_n(n, a, lda, ipiv, work, lwork, info, dc_access); \
  2324. break; \
  2325. } \
  2326. } \
  2327. MKL_DC_CONVERT_INT(work[0], MKL_DC_MAX(1, n)); \
  2328. } while (0)
  2329. static __inline void mkl_dc_getri(MKL_INT n, mkl_dc_type* a, MKL_INT lda, const MKL_INT* ipiv,
  2330. mkl_dc_type* work, MKL_INT lwork, MKL_INT* info)
  2331. {
  2332. mkl_dc_getri_generic(n, a, lda, ipiv, work, lwork, info, MKL_DC_MN);
  2333. }
  2334. #ifndef MKL_DIRECT_CALL_LAPACKE_DISABLE
  2335. static __inline lapack_int mkl_dc_lapacke_getri_convert(int matrix_layout, lapack_int n, mkl_dc_type* a,
  2336. lapack_int lda, const lapack_int* ipiv,
  2337. mkl_dc_type* work, lapack_int lwork)
  2338. {
  2339. lapack_int info = 0;
  2340. if (MKL_DC_GETRI_CHECKSIZE(n)) {
  2341. mkl_dc_type work_small[10];
  2342. lapack_int lwork_small = 10;
  2343. if (!work) {
  2344. work = work_small;
  2345. lwork = lwork_small;
  2346. }
  2347. if (matrix_layout == LAPACK_ROW_MAJOR) {
  2348. mkl_dc_getri_generic(n, a, lda, ipiv, work, lwork, &info, MKL_DC_MT);
  2349. } else {
  2350. mkl_dc_getri_generic(n, a, lda, ipiv, work, lwork, &info, MKL_DC_MN);
  2351. }
  2352. return info;
  2353. }
  2354. if (work) {
  2355. return MKL_DC_CONCAT3(LAPACKE_, MKL_DC_PREC_LETTER, getri_work)(matrix_layout, n, a, lda, ipiv, work, lwork);
  2356. }
  2357. return MKL_DC_CONCAT3(LAPACKE_, MKL_DC_PREC_LETTER, getri)(matrix_layout, n, a, lda, ipiv);
  2358. }
  2359. #endif
  2360. /* ?GEQRF */
  2361. #if defined(MKL_SINGLE) || defined(MKL_DOUBLE)
  2362. #define mkl_dc_geqrf_hh(n, alpha, x, ldx, tau, dc_access) \
  2363. do { \
  2364. MKL_INT ii; \
  2365. MKL_INT n1 = n - 1; \
  2366. mkl_dc_real_type xnorm; \
  2367. mkl_dc_real_type ax; \
  2368. if (n <= 1) { \
  2369. MKL_DC_SET_ZERO(tau); \
  2370. break; \
  2371. } \
  2372. MKL_DC_NRM2_VEC(xnorm, n1, x, ldx, dc_access); \
  2373. if (MKL_DC_IS_R_ZERO(xnorm)) { \
  2374. MKL_DC_SET_ZERO(tau); \
  2375. } else { \
  2376. MKL_INT knt = 0; \
  2377. mkl_dc_type beta; \
  2378. mkl_dc_real_type safmin; \
  2379. mkl_dc_real_type ab; \
  2380. MKL_DC_R_PY2(ax, alpha, xnorm); \
  2381. MKL_DC_R_SIGN(beta, ax, alpha); \
  2382. MKL_DC_NEG(beta, beta); \
  2383. if (!MKL_DC_UNSAFE) { \
  2384. safmin = MKL_DC_XLAMCH("s") / MKL_DC_XLAMCH("e"); \
  2385. if (MKL_DC_ABS(beta) < safmin) { \
  2386. mkl_dc_real_type rsafmn; \
  2387. MKL_DC_INV(rsafmn, safmin); \
  2388. do { \
  2389. knt++; \
  2390. for (ii = 0; ii < n - 1; ++ii) { \
  2391. MKL_DC_MUL(dc_access(x, ldx, ii, 0), dc_access(x, ldx, ii, 0), rsafmn); \
  2392. } \
  2393. MKL_DC_MUL(beta, beta, rsafmn); \
  2394. MKL_DC_MUL(alpha, alpha, rsafmn); \
  2395. } while (MKL_DC_ABS(beta) < safmin); \
  2396. MKL_DC_NRM2_VEC(xnorm, n1, x, ldx, dc_access); \
  2397. MKL_DC_R_PY2(ax, alpha, xnorm); \
  2398. MKL_DC_R_SIGN(beta, ax, alpha); \
  2399. MKL_DC_R_NEG(beta, beta); \
  2400. } \
  2401. } \
  2402. MKL_DC_SUB(tau, beta, alpha); \
  2403. MKL_DC_DIV(tau, tau, beta); \
  2404. MKL_DC_SUB(ab, alpha, beta); \
  2405. for (ii = 0; ii < n - 1; ++ii) { \
  2406. MKL_DC_DIV(dc_access(x, ldx, ii, 0), dc_access(x, ldx, ii, 0), ab); \
  2407. } \
  2408. if (!MKL_DC_UNSAFE) { \
  2409. for (ii = 0; ii < knt; ++ii) { \
  2410. MKL_DC_R_MUL(beta, beta, safmin); \
  2411. } \
  2412. } \
  2413. alpha = beta; \
  2414. } \
  2415. } while (0)
  2416. #endif
  2417. #if defined(MKL_COMPLEX) || defined(MKL_COMPLEX16)
  2418. #define mkl_dc_geqrf_hh(n, alpha, x, ldx, tau, dc_access) \
  2419. do { \
  2420. MKL_INT ii; \
  2421. MKL_INT n1 = n - 1; \
  2422. mkl_dc_real_type xnorm; \
  2423. mkl_dc_real_type ax; \
  2424. mkl_dc_real_type alphr; \
  2425. mkl_dc_real_type alphi; \
  2426. mkl_dc_real_type beta; \
  2427. mkl_dc_real_type rzero; \
  2428. mkl_dc_type temp; \
  2429. mkl_dc_real_type re_temp; \
  2430. mkl_dc_real_type im_temp; \
  2431. MKL_DC_SET_R_ZERO(rzero); \
  2432. if (n <= 0) { \
  2433. MKL_DC_SET_ZERO(tau); \
  2434. break; \
  2435. } \
  2436. MKL_DC_NRM2_VEC(xnorm, n1, x, ldx, dc_access); \
  2437. MKL_DC_REIM(alphr, alphi, alpha); \
  2438. if (MKL_DC_IS_R_ZERO(xnorm) && MKL_DC_IS_R_ZERO(alphi)) { \
  2439. MKL_DC_SET_ZERO(tau); \
  2440. } else { \
  2441. MKL_INT knt = 0; \
  2442. mkl_dc_real_type safmin; \
  2443. MKL_DC_R_PY3(ax, alphr, alphi, xnorm); \
  2444. MKL_DC_R_SIGN(beta, ax, alphr); \
  2445. MKL_DC_R_NEG(beta, beta); \
  2446. if (!MKL_DC_UNSAFE) { \
  2447. safmin = MKL_DC_XLAMCH("s") / MKL_DC_XLAMCH("e"); \
  2448. if (MKL_DC_ABS(beta) < safmin) { \
  2449. mkl_dc_real_type rsafmn; \
  2450. mkl_dc_type c_rsafmn; \
  2451. MKL_DC_R_INV(rsafmn, safmin); \
  2452. MKL_DC_CMPLX(c_rsafmn, rsafmn, rzero); \
  2453. do { \
  2454. knt++; \
  2455. for (ii = 0; ii < n - 1; ++ii) { \
  2456. MKL_DC_MUL(dc_access(x, ldx, ii, 0), dc_access(x, ldx, ii, 0), c_rsafmn); \
  2457. } \
  2458. MKL_DC_R_MUL(beta, beta, rsafmn); \
  2459. MKL_DC_R_MUL(alphi, alphi, rsafmn); \
  2460. MKL_DC_R_MUL(alphr, alphr, rsafmn); \
  2461. } while (MKL_DC_ABS(beta) < safmin); \
  2462. MKL_DC_NRM2_VEC(xnorm, n1, x, ldx, dc_access); \
  2463. MKL_DC_CMPLX(alpha, alphr, alphi); \
  2464. MKL_DC_R_PY3(ax, alphr, alphi, xnorm); \
  2465. MKL_DC_R_SIGN(beta, ax, alphr); \
  2466. MKL_DC_R_NEG(beta, beta); \
  2467. } \
  2468. } \
  2469. MKL_DC_R_SUB(re_temp, beta, alphr); \
  2470. MKL_DC_R_DIV(re_temp, re_temp, beta); \
  2471. MKL_DC_R_NEG(im_temp, alphi); \
  2472. MKL_DC_R_DIV(im_temp, im_temp, beta); \
  2473. MKL_DC_CMPLX(tau, re_temp, im_temp); \
  2474. MKL_DC_CMPLX(temp, beta, rzero); \
  2475. MKL_DC_SUB(temp, alpha, temp); \
  2476. MKL_DC_INV(alpha, temp); \
  2477. for (ii = 0; ii < n - 1; ++ii) { \
  2478. MKL_DC_MUL(dc_access(x, ldx, ii, 0), dc_access(x, ldx, ii, 0), alpha); \
  2479. } \
  2480. if (!MKL_DC_UNSAFE) { \
  2481. for (ii = 0; ii < knt; ++ii) { \
  2482. MKL_DC_R_MUL(beta, beta, safmin); \
  2483. } \
  2484. } \
  2485. MKL_DC_CMPLX(alpha, beta, rzero); \
  2486. } \
  2487. } while (0)
  2488. #endif
  2489. #define mkl_dc_geqrf_apply(m, n, v, ldv, tau, c, ldc, work, dc_access) \
  2490. do { \
  2491. MKL_INT ii; \
  2492. MKL_INT jj; \
  2493. mkl_dc_type mtau; \
  2494. MKL_DC_NEG(mtau, tau); \
  2495. for (ii = 0; ii < n; ++ii) { \
  2496. MKL_DC_SET_ZERO(MKL_DC_ACCESS1D(work, ii)); \
  2497. for (jj = 0; jj < m; ++jj) { \
  2498. mkl_dc_type cji; \
  2499. MKL_DC_CONJ(cji, dc_access(c, ldc, jj, ii)); \
  2500. MKL_DC_MUL_ADD(MKL_DC_ACCESS1D(work, ii), cji, dc_access(v, ldv, jj, 0), MKL_DC_ACCESS1D(work, ii)); \
  2501. } \
  2502. } \
  2503. for (jj = 0; jj < n; ++jj) { \
  2504. mkl_dc_type temp; \
  2505. mkl_dc_type workj; \
  2506. MKL_DC_CONJ(workj, MKL_DC_ACCESS1D(work, jj)); \
  2507. MKL_DC_MUL(temp, mtau, workj); \
  2508. for (ii = 0; ii < m; ++ii) { \
  2509. MKL_DC_MUL_ADD(dc_access(c, ldc, ii, jj), dc_access(v, ldv, ii, 0), temp, dc_access(c, ldc, ii, jj)); \
  2510. } \
  2511. } \
  2512. } while (0)
  2513. #define mkl_dc_geqrf_mxn(m, n, a, lda, tau, work, lwork, info, dc_access) \
  2514. do { \
  2515. MKL_INT i; \
  2516. MKL_INT k = MKL_DC_MIN(m, n); \
  2517. if (lwork == -1) { \
  2518. MKL_DC_CONVERT_INT(MKL_DC_ACCESS1D(work, 0), MKL_DC_MAX(1, n)); \
  2519. } \
  2520. for (i = 0; i < k; ++i) { \
  2521. MKL_INT i1 = MKL_DC_MIN(i + 1, m - 1); \
  2522. MKL_INT mi = m - i; \
  2523. MKL_INT ni1 = n - i - 1; \
  2524. mkl_dc_geqrf_hh((m - i), (dc_access(a, lda, i, i)), (&dc_access(a, lda, i1, i)), lda, MKL_DC_ACCESS1D(tau, i), dc_access); \
  2525. if (i + 1 < n) { \
  2526. mkl_dc_type aii = dc_access(a, lda, i, i); \
  2527. mkl_dc_type taui; \
  2528. MKL_DC_SET_ONE(dc_access(a, lda, i, i)); \
  2529. MKL_DC_CONJ(taui, MKL_DC_ACCESS1D(tau, i)); \
  2530. mkl_dc_geqrf_apply(mi, ni1, (&dc_access(a, lda, i, i)), lda, taui, (&dc_access(a, lda, i, i + 1)), lda, work, dc_access); \
  2531. dc_access(a, lda, i, i) = aii; \
  2532. } \
  2533. } \
  2534. MKL_DC_CONVERT_INT(MKL_DC_ACCESS1D(work, 0), MKL_DC_MAX(1, n)); \
  2535. } while (0)
  2536. #define mkl_dc_geqrf_generic(m, n, a, lda, tau, work, lwork, info, dc_access) \
  2537. do { \
  2538. *info = 0; \
  2539. if (lwork != -1) { \
  2540. if (m == n && m <= 5) { \
  2541. switch (n) { \
  2542. case 0: \
  2543. break; \
  2544. case 1: \
  2545. mkl_dc_geqrf_mxn(1, 1, a, lda, tau, work, lwork, info, dc_access); \
  2546. break; \
  2547. case 2: \
  2548. mkl_dc_geqrf_mxn(2, 2, a, lda, tau, work, lwork, info, dc_access); \
  2549. break; \
  2550. case 3: \
  2551. mkl_dc_geqrf_mxn(3, 3, a, lda, tau, work, lwork, info, dc_access); \
  2552. break; \
  2553. case 4: \
  2554. mkl_dc_geqrf_mxn(4, 4, a, lda, tau, work, lwork, info, dc_access); \
  2555. break; \
  2556. case 5: \
  2557. mkl_dc_geqrf_mxn(5, 5, a, lda, tau, work, lwork, info, dc_access); \
  2558. break; \
  2559. } \
  2560. } else { \
  2561. mkl_dc_geqrf_mxn(m, n, a, lda, tau, work, lwork, info, dc_access); \
  2562. } \
  2563. } else { \
  2564. MKL_DC_CONVERT_INT(MKL_DC_ACCESS1D(work, 0), MKL_DC_MAX(1, n)); \
  2565. } \
  2566. } while (0)
  2567. static __inline void mkl_dc_geqrf(MKL_INT m, MKL_INT n, mkl_dc_type* a, MKL_INT lda, mkl_dc_type* tau,
  2568. mkl_dc_type* work, MKL_INT lwork, MKL_INT* info)
  2569. {
  2570. mkl_dc_geqrf_generic(m, n, a, lda, tau, work, lwork, info, MKL_DC_MN);
  2571. }
  2572. #ifndef MKL_DIRECT_CALL_LAPACKE_DISABLE
  2573. static __inline lapack_int mkl_dc_lapacke_geqrf_convert(int matrix_layout, lapack_int m, lapack_int n,
  2574. mkl_dc_type* a, lapack_int lda,
  2575. mkl_dc_type* tau, mkl_dc_type* work, lapack_int lwork)
  2576. {
  2577. lapack_int info = 0;
  2578. if (MKL_DC_GEQRF_CHECKSIZE(m, n)) {
  2579. mkl_dc_type work_small[10];
  2580. lapack_int lwork_small = 10;
  2581. if (!work) {
  2582. work = work_small;
  2583. lwork = lwork_small;
  2584. }
  2585. if (matrix_layout == LAPACK_ROW_MAJOR) {
  2586. mkl_dc_geqrf_generic(m, n, a, lda, tau, work, lwork, &info, MKL_DC_MT);
  2587. } else {
  2588. mkl_dc_geqrf_generic(m, n, a, lda, tau, work, lwork, &info, MKL_DC_MN);
  2589. }
  2590. return info;
  2591. }
  2592. if (work) {
  2593. return MKL_DC_CONCAT3(LAPACKE_, MKL_DC_PREC_LETTER, geqrf_work)(matrix_layout, m, n, a, lda, tau, work, lwork);
  2594. }
  2595. return MKL_DC_CONCAT3(LAPACKE_, MKL_DC_PREC_LETTER, geqrf)(matrix_layout, m, n, a, lda, tau);
  2596. }
  2597. #endif
  2598. /* ?POTRF */
  2599. #define mkl_dc_potrf_n(uplo, n, a, lda, info, dc_access) \
  2600. do { \
  2601. MKL_INT i, j, ii, jj; \
  2602. if (MKL_DC_MisU(uplo)) { \
  2603. for (j = 0; j < n; ++j) { \
  2604. mkl_dc_type ajj = dc_access(a, lda, j, j); \
  2605. MKL_DC_ZERO_IMAG(ajj); \
  2606. for (i = 0; i < j; ++i) { \
  2607. MKL_DC_SUB_MUL_CONJ(ajj, ajj, dc_access(a, lda, i, j), dc_access(a, lda, i, j)); \
  2608. } \
  2609. if (MKL_DC_NON_POS(ajj)) { \
  2610. dc_access(a, lda, j, j) = ajj; \
  2611. *info = j + 1; \
  2612. break; \
  2613. } \
  2614. MKL_DC_SQRT(ajj, ajj); \
  2615. dc_access(a, lda, j, j) = ajj; \
  2616. MKL_DC_INV(ajj, ajj); \
  2617. if (j < n - 1) { \
  2618. for (jj = j + 1; jj < n; ++jj) { \
  2619. mkl_dc_type aj_jj = dc_access(a, lda, j, jj); \
  2620. for (ii = 0; ii < j; ++ii) { \
  2621. MKL_DC_SUB_MUL_CONJ(aj_jj, aj_jj, \
  2622. dc_access(a, lda, ii, jj), dc_access(a, lda, ii, j)); \
  2623. } \
  2624. dc_access(a, lda, j, jj) = aj_jj; \
  2625. } \
  2626. for (i = j + 1; i < n; ++i) { \
  2627. MKL_DC_MUL(dc_access(a, lda, j, i), dc_access(a, lda, j, i), ajj); \
  2628. } \
  2629. } \
  2630. } \
  2631. } else { \
  2632. for (j = 0; j < n; ++j) { \
  2633. mkl_dc_type ajj = dc_access(a, lda, j, j); \
  2634. MKL_DC_ZERO_IMAG(ajj); \
  2635. for (i = 0; i < j; ++i) { \
  2636. MKL_DC_SUB_MUL_CONJ(ajj, ajj, dc_access(a, lda, j, i), dc_access(a, lda, j, i)); \
  2637. } \
  2638. if (MKL_DC_NON_POS(ajj)) { \
  2639. dc_access(a, lda, j, j) = ajj; \
  2640. *info = j + 1; \
  2641. break; \
  2642. } \
  2643. MKL_DC_SQRT(ajj, ajj); \
  2644. dc_access(a, lda, j, j) = ajj; \
  2645. MKL_DC_INV(ajj, ajj); \
  2646. if (j < n - 1) { \
  2647. for (ii = 0; ii < j; ++ii) { \
  2648. for (jj = j + 1; jj < n; ++jj) { \
  2649. MKL_DC_SUB_MUL_CONJ(dc_access(a, lda, jj, j), dc_access(a, lda, jj, j), \
  2650. dc_access(a, lda, jj, ii), dc_access(a, lda, j, ii)); \
  2651. } \
  2652. } \
  2653. for (i = j + 1; i < n; ++i) { \
  2654. MKL_DC_MUL(dc_access(a, lda, i, j), dc_access(a, lda, i, j), ajj); \
  2655. } \
  2656. } \
  2657. } \
  2658. } \
  2659. } while (0) \
  2660. #define mkl_dc_potrf_1(uplo, n, a, lda, info, dc_access) \
  2661. do { \
  2662. mkl_dc_type a00 = dc_access(a, lda, 0, 0); \
  2663. MKL_DC_ZERO_IMAG(a00); \
  2664. if (MKL_DC_NON_POS(a00)) { \
  2665. dc_access(a, lda, 0, 0) = a00; \
  2666. *info = 1; \
  2667. break; \
  2668. } \
  2669. MKL_DC_SQRT(a00, a00); \
  2670. dc_access(a, lda, 0, 0) = a00; \
  2671. } while (0) \
  2672. #define mkl_dc_potrf_2(uplo, n, a, lda, info, dc_access) \
  2673. do { \
  2674. mkl_dc_type a00 = dc_access(a, lda, 0, 0); \
  2675. mkl_dc_type a11 = dc_access(a, lda, 1, 1); \
  2676. MKL_DC_ZERO_IMAG(a00); \
  2677. MKL_DC_ZERO_IMAG(a11); \
  2678. if (MKL_DC_MisU(uplo)) { \
  2679. mkl_dc_type a01 = dc_access(a, lda, 0, 1); \
  2680. if (MKL_DC_NON_POS(a00)) { \
  2681. dc_access(a, lda, 0, 0) = a00; \
  2682. *info = 1; \
  2683. break; \
  2684. } \
  2685. MKL_DC_SQRT(a00, a00); \
  2686. dc_access(a, lda, 0, 0) = a00; \
  2687. MKL_DC_INV(a00, a00); \
  2688. MKL_DC_MUL(a01, a01, a00); \
  2689. dc_access(a, lda, 0, 1) = a01; \
  2690. MKL_DC_SUB_MUL_CONJ(a11, a11, a01, a01); \
  2691. if (MKL_DC_NON_POS(a11)) { \
  2692. dc_access(a, lda, 1, 1) = a11; \
  2693. *info = 2; \
  2694. break; \
  2695. } \
  2696. MKL_DC_SQRT(a11, a11); \
  2697. dc_access(a, lda, 1, 1) = a11; \
  2698. } else { \
  2699. mkl_dc_type a10 = dc_access(a, lda, 1, 0); \
  2700. if (MKL_DC_NON_POS(a00)) { \
  2701. dc_access(a, lda, 0, 0) = a00; \
  2702. *info = 1; \
  2703. break; \
  2704. } \
  2705. MKL_DC_SQRT(a00, a00); \
  2706. dc_access(a, lda, 0, 0) = a00; \
  2707. MKL_DC_INV(a00, a00); \
  2708. MKL_DC_MUL(a10, a10, a00); \
  2709. dc_access(a, lda, 1, 0) = a10; \
  2710. MKL_DC_SUB_MUL_CONJ(a11, a11, a10, a10); \
  2711. if (MKL_DC_NON_POS(a11)) { \
  2712. dc_access(a, lda, 1, 1) = a11; \
  2713. *info = 2; \
  2714. break; \
  2715. } \
  2716. MKL_DC_SQRT(a11, a11); \
  2717. dc_access(a, lda, 1, 1) = a11; \
  2718. } \
  2719. } while (0)
  2720. #define mkl_dc_potrf_3(uplo, n, a, lda, info, dc_access) \
  2721. do { \
  2722. mkl_dc_type a00 = dc_access(a, lda, 0, 0); \
  2723. mkl_dc_type a11 = dc_access(a, lda, 1, 1); \
  2724. mkl_dc_type a22 = dc_access(a, lda, 2, 2); \
  2725. MKL_DC_ZERO_IMAG(a00); \
  2726. MKL_DC_ZERO_IMAG(a11); \
  2727. MKL_DC_ZERO_IMAG(a22); \
  2728. if (MKL_DC_MisU(uplo)) { \
  2729. mkl_dc_type a01 = dc_access(a, lda, 0, 1); \
  2730. mkl_dc_type a12 = dc_access(a, lda, 1, 2); \
  2731. mkl_dc_type a02 = dc_access(a, lda, 0, 2); \
  2732. if (MKL_DC_NON_POS(a00)) { \
  2733. dc_access(a, lda, 0, 0) = a00; \
  2734. *info = 1; \
  2735. break; \
  2736. } \
  2737. MKL_DC_SQRT(a00, a00); \
  2738. dc_access(a, lda, 0, 0) = a00; \
  2739. MKL_DC_INV(a00, a00); \
  2740. MKL_DC_MUL(a01, a01, a00); \
  2741. dc_access(a, lda, 0, 1) = a01; \
  2742. MKL_DC_MUL(a02, a02, a00); \
  2743. dc_access(a, lda, 0, 2) = a02; \
  2744. MKL_DC_SUB_MUL_CONJ(a11, a11, a01, a01); \
  2745. if (MKL_DC_NON_POS(a11)) { \
  2746. dc_access(a, lda, 1, 1) = a11; \
  2747. *info = 2; \
  2748. break; \
  2749. } \
  2750. MKL_DC_SQRT(a11, a11); \
  2751. dc_access(a, lda, 1, 1) = a11; \
  2752. MKL_DC_SUB_MUL_CONJ(a12, a12, a02, a01); \
  2753. MKL_DC_INV(a11, a11); \
  2754. MKL_DC_MUL(a12, a12, a11); \
  2755. dc_access(a, lda, 1, 2) = a12; \
  2756. MKL_DC_SUB_MUL_CONJ(a22, a22, a02, a02); \
  2757. MKL_DC_SUB_MUL_CONJ(a22, a22, a12, a12); \
  2758. if (MKL_DC_NON_POS(a22)) { \
  2759. dc_access(a, lda, 2, 2) = a22; \
  2760. *info = 3; \
  2761. break; \
  2762. } \
  2763. MKL_DC_SQRT(a22, a22); \
  2764. dc_access(a, lda, 2, 2) = a22; \
  2765. } else { \
  2766. mkl_dc_type a10 = dc_access(a, lda, 1, 0); \
  2767. mkl_dc_type a20 = dc_access(a, lda, 2, 0); \
  2768. mkl_dc_type a21 = dc_access(a, lda, 2, 1); \
  2769. if (MKL_DC_NON_POS(a00)) { \
  2770. dc_access(a, lda, 0, 0) = a00; \
  2771. *info = 1; \
  2772. break; \
  2773. } \
  2774. MKL_DC_SQRT(a00, a00); \
  2775. dc_access(a, lda, 0, 0) = a00; \
  2776. MKL_DC_INV(a00, a00); \
  2777. MKL_DC_MUL(a10, a10, a00); \
  2778. dc_access(a, lda, 1, 0) = a10; \
  2779. MKL_DC_MUL(a20, a20, a00); \
  2780. dc_access(a, lda, 2, 0) = a20; \
  2781. MKL_DC_SUB_MUL_CONJ(a11, a11, a10, a10); \
  2782. if (MKL_DC_NON_POS(a11)) { \
  2783. dc_access(a, lda, 1, 1) = a11; \
  2784. *info = 2; \
  2785. break; \
  2786. } \
  2787. MKL_DC_SQRT(a11, a11); \
  2788. dc_access(a, lda, 1, 1) = a11; \
  2789. MKL_DC_SUB_MUL_CONJ(a21, a21, a20, a10); \
  2790. MKL_DC_INV(a11, a11); \
  2791. MKL_DC_MUL(a21, a21, a11); \
  2792. dc_access(a, lda, 2, 1) = a21; \
  2793. MKL_DC_SUB_MUL_CONJ(a22, a22, a20, a20); \
  2794. MKL_DC_SUB_MUL_CONJ(a22, a22, a21, a21); \
  2795. if (MKL_DC_NON_POS(a22)) { \
  2796. dc_access(a, lda, 2, 2) = a22; \
  2797. *info = 3; \
  2798. break; \
  2799. } \
  2800. MKL_DC_SQRT(a22, a22); \
  2801. dc_access(a, lda, 2, 2) = a22; \
  2802. } \
  2803. } while (0)
  2804. #define mkl_dc_potrf_4(uplo, n, a, lda, info, dc_access) \
  2805. do { \
  2806. mkl_dc_type a00 = dc_access(a, lda, 0, 0); \
  2807. mkl_dc_type a11 = dc_access(a, lda, 1, 1); \
  2808. mkl_dc_type a22 = dc_access(a, lda, 2, 2); \
  2809. mkl_dc_type a33 = dc_access(a, lda, 3, 3); \
  2810. MKL_DC_ZERO_IMAG(a00); \
  2811. MKL_DC_ZERO_IMAG(a11); \
  2812. MKL_DC_ZERO_IMAG(a22); \
  2813. MKL_DC_ZERO_IMAG(a33); \
  2814. if (MKL_DC_MisU(uplo)) { \
  2815. mkl_dc_type a01 = dc_access(a, lda, 0, 1); \
  2816. mkl_dc_type a02 = dc_access(a, lda, 0, 2); \
  2817. mkl_dc_type a03 = dc_access(a, lda, 0, 3); \
  2818. mkl_dc_type a12 = dc_access(a, lda, 1, 2); \
  2819. mkl_dc_type a13 = dc_access(a, lda, 1, 3); \
  2820. mkl_dc_type a23 = dc_access(a, lda, 2, 3); \
  2821. if (MKL_DC_NON_POS(a00)) { \
  2822. dc_access(a, lda, 0, 0) = a00; \
  2823. *info = 1; \
  2824. break; \
  2825. } \
  2826. MKL_DC_SQRT(a00, a00); \
  2827. dc_access(a, lda, 0, 0) = a00; \
  2828. MKL_DC_INV(a00, a00); \
  2829. MKL_DC_MUL(a01, a01, a00); \
  2830. dc_access(a, lda, 0, 1) = a01; \
  2831. MKL_DC_MUL(a02, a02, a00); \
  2832. dc_access(a, lda, 0, 2) = a02; \
  2833. MKL_DC_MUL(a03, a03, a00); \
  2834. dc_access(a, lda, 0, 3) = a03; \
  2835. MKL_DC_SUB_MUL_CONJ(a11, a11, a01, a01); \
  2836. if (MKL_DC_NON_POS(a11)) { \
  2837. dc_access(a, lda, 1, 1) = a11; \
  2838. *info = 2; \
  2839. break; \
  2840. } \
  2841. MKL_DC_SQRT(a11, a11); \
  2842. dc_access(a, lda, 1, 1) = a11; \
  2843. MKL_DC_SUB_MUL_CONJ(a12, a12, a02, a01); \
  2844. MKL_DC_SUB_MUL_CONJ(a13, a13, a03, a01); \
  2845. MKL_DC_INV(a11, a11); \
  2846. MKL_DC_MUL(a12, a12, a11); \
  2847. dc_access(a, lda, 1, 2) = a12; \
  2848. MKL_DC_MUL(a13, a13, a11); \
  2849. dc_access(a, lda, 1, 3) = a13; \
  2850. MKL_DC_SUB_MUL_CONJ(a22, a22, a02, a02); \
  2851. MKL_DC_SUB_MUL_CONJ(a22, a22, a12, a12); \
  2852. if (MKL_DC_NON_POS(a22)) { \
  2853. dc_access(a, lda, 2, 2) = a22; \
  2854. *info = 3; \
  2855. break; \
  2856. } \
  2857. MKL_DC_SQRT(a22, a22); \
  2858. dc_access(a, lda, 2, 2) = a22; \
  2859. MKL_DC_SUB_MUL_CONJ(a23, a23, a03, a02); \
  2860. MKL_DC_SUB_MUL_CONJ(a23, a23, a13, a12); \
  2861. MKL_DC_INV(a22, a22); \
  2862. MKL_DC_MUL(a23, a23, a22); \
  2863. dc_access(a, lda, 2, 3) = a23; \
  2864. MKL_DC_SUB_MUL_CONJ(a33, a33, a03, a03); \
  2865. MKL_DC_SUB_MUL_CONJ(a33, a33, a13, a13); \
  2866. MKL_DC_SUB_MUL_CONJ(a33, a33, a23, a23); \
  2867. if (MKL_DC_NON_POS(a33)) { \
  2868. dc_access(a, lda, 3, 3) = a33; \
  2869. *info = 4; \
  2870. break; \
  2871. } \
  2872. MKL_DC_SQRT(a33, a33); \
  2873. dc_access(a, lda, 3, 3) = a33; \
  2874. } else { \
  2875. mkl_dc_type a10 = dc_access(a, lda, 1, 0); \
  2876. mkl_dc_type a20 = dc_access(a, lda, 2, 0); \
  2877. mkl_dc_type a30 = dc_access(a, lda, 3, 0); \
  2878. mkl_dc_type a21 = dc_access(a, lda, 2, 1); \
  2879. mkl_dc_type a31 = dc_access(a, lda, 3, 1); \
  2880. mkl_dc_type a32 = dc_access(a, lda, 3, 2); \
  2881. if (MKL_DC_NON_POS(a00)) { \
  2882. dc_access(a, lda, 0, 0) = a00; \
  2883. *info = 1; \
  2884. break; \
  2885. } \
  2886. MKL_DC_SQRT(a00, a00); \
  2887. dc_access(a, lda, 0, 0) = a00; \
  2888. MKL_DC_INV(a00, a00); \
  2889. MKL_DC_MUL(a10, a10, a00); \
  2890. dc_access(a, lda, 1, 0) = a10; \
  2891. MKL_DC_MUL(a20, a20, a00); \
  2892. dc_access(a, lda, 2, 0) = a20; \
  2893. MKL_DC_MUL(a30, a30, a00); \
  2894. dc_access(a, lda, 3, 0) = a30; \
  2895. MKL_DC_SUB_MUL_CONJ(a11, a11, a10, a10); \
  2896. if (MKL_DC_NON_POS(a11)) { \
  2897. dc_access(a, lda, 1, 1) = a11; \
  2898. *info = 2; \
  2899. break; \
  2900. } \
  2901. MKL_DC_SQRT(a11, a11); \
  2902. dc_access(a, lda, 1, 1) = a11; \
  2903. MKL_DC_SUB_MUL_CONJ(a21, a21, a20, a10); \
  2904. MKL_DC_SUB_MUL_CONJ(a31, a31, a30, a10); \
  2905. MKL_DC_INV(a11, a11); \
  2906. MKL_DC_MUL(a21, a21, a11); \
  2907. dc_access(a, lda, 2, 1) = a21; \
  2908. MKL_DC_MUL(a31, a31, a11); \
  2909. dc_access(a, lda, 3, 1) = a31; \
  2910. MKL_DC_SUB_MUL_CONJ(a22, a22, a20, a20); \
  2911. MKL_DC_SUB_MUL_CONJ(a22, a22, a21, a21); \
  2912. if (MKL_DC_NON_POS(a22)) { \
  2913. dc_access(a, lda, 2, 2) = a22; \
  2914. *info = 3; \
  2915. break; \
  2916. } \
  2917. MKL_DC_SQRT(a22, a22); \
  2918. dc_access(a, lda, 2, 2) = a22; \
  2919. MKL_DC_SUB_MUL_CONJ(a32, a32, a30, a20); \
  2920. MKL_DC_SUB_MUL_CONJ(a32, a32, a31, a21); \
  2921. MKL_DC_INV(a22, a22); \
  2922. MKL_DC_MUL(a32, a32, a22); \
  2923. dc_access(a, lda, 3, 2) = a32; \
  2924. MKL_DC_SUB_MUL_CONJ(a33, a33, a30, a30); \
  2925. MKL_DC_SUB_MUL_CONJ(a33, a33, a31, a31); \
  2926. MKL_DC_SUB_MUL_CONJ(a33, a33, a32, a32); \
  2927. if (MKL_DC_NON_POS(a33)) { \
  2928. dc_access(a, lda, 3, 3) = a33; \
  2929. *info = 4; \
  2930. break; \
  2931. } \
  2932. MKL_DC_SQRT(a33, a33); \
  2933. dc_access(a, lda, 3, 3) = a33; \
  2934. } \
  2935. } while (0)
  2936. #define mkl_dc_potrf_5(uplo, n, a, lda, info, dc_access) \
  2937. { \
  2938. mkl_dc_type a00 = dc_access(a, lda, 0, 0); \
  2939. mkl_dc_type a11 = dc_access(a, lda, 1, 1); \
  2940. mkl_dc_type a22 = dc_access(a, lda, 2, 2); \
  2941. mkl_dc_type a33 = dc_access(a, lda, 3, 3); \
  2942. mkl_dc_type a44 = dc_access(a, lda, 4, 4); \
  2943. MKL_DC_ZERO_IMAG(a00); \
  2944. MKL_DC_ZERO_IMAG(a11); \
  2945. MKL_DC_ZERO_IMAG(a22); \
  2946. MKL_DC_ZERO_IMAG(a33); \
  2947. MKL_DC_ZERO_IMAG(a44); \
  2948. if (MKL_DC_MisU(uplo)) { \
  2949. mkl_dc_type a01 = dc_access(a, lda, 0, 1); \
  2950. mkl_dc_type a02 = dc_access(a, lda, 0, 2); \
  2951. mkl_dc_type a12 = dc_access(a, lda, 1, 2); \
  2952. mkl_dc_type a03 = dc_access(a, lda, 0, 3); \
  2953. mkl_dc_type a13 = dc_access(a, lda, 1, 3); \
  2954. mkl_dc_type a23 = dc_access(a, lda, 2, 3); \
  2955. mkl_dc_type a04 = dc_access(a, lda, 0, 4); \
  2956. mkl_dc_type a14 = dc_access(a, lda, 1, 4); \
  2957. mkl_dc_type a24 = dc_access(a, lda, 2, 4); \
  2958. mkl_dc_type a34 = dc_access(a, lda, 3, 4); \
  2959. if (MKL_DC_NON_POS(a00)) { \
  2960. dc_access(a, lda, 0, 0) = a00; \
  2961. *info = 1; \
  2962. break; \
  2963. } \
  2964. MKL_DC_SQRT(a00, a00); \
  2965. dc_access(a, lda, 0, 0) = a00; \
  2966. MKL_DC_INV(a00, a00); \
  2967. MKL_DC_MUL(a01, a01, a00); \
  2968. MKL_DC_MUL(a02, a02, a00); \
  2969. MKL_DC_MUL(a03, a03, a00); \
  2970. MKL_DC_MUL(a04, a04, a00); \
  2971. dc_access(a, lda, 0, 1) = a01; \
  2972. dc_access(a, lda, 0, 2) = a02; \
  2973. dc_access(a, lda, 0, 3) = a03; \
  2974. dc_access(a, lda, 0, 4) = a04; \
  2975. MKL_DC_SUB_MUL_CONJ(a11, a11, a01, a01); \
  2976. if (MKL_DC_NON_POS(a11)) { \
  2977. dc_access(a, lda, 1, 1) = a11; \
  2978. *info = 2; \
  2979. break; \
  2980. } \
  2981. MKL_DC_SQRT(a11, a11); \
  2982. dc_access(a, lda, 1, 1) = a11; \
  2983. MKL_DC_SUB_MUL_CONJ(a12, a12, a02, a01); \
  2984. MKL_DC_SUB_MUL_CONJ(a13, a13, a03, a01); \
  2985. MKL_DC_SUB_MUL_CONJ(a14, a14, a04, a01); \
  2986. MKL_DC_INV(a11, a11); \
  2987. MKL_DC_MUL(a12, a12, a11); \
  2988. MKL_DC_MUL(a13, a13, a11); \
  2989. MKL_DC_MUL(a14, a14, a11); \
  2990. dc_access(a, lda, 1, 2) = a12; \
  2991. dc_access(a, lda, 1, 3) = a13; \
  2992. dc_access(a, lda, 1, 4) = a14; \
  2993. MKL_DC_SUB_MUL_CONJ(a22, a22, a02, a02); \
  2994. MKL_DC_SUB_MUL_CONJ(a22, a22, a12, a12); \
  2995. if (MKL_DC_NON_POS(a22)) { \
  2996. dc_access(a, lda, 2, 2) = a22; \
  2997. *info = 3; \
  2998. break; \
  2999. } \
  3000. MKL_DC_SQRT(a22, a22); \
  3001. dc_access(a, lda, 2, 2) = a22; \
  3002. MKL_DC_SUB_MUL_CONJ(a23, a23, a03, a02); \
  3003. MKL_DC_SUB_MUL_CONJ(a23, a23, a13, a12); \
  3004. MKL_DC_SUB_MUL_CONJ(a24, a24, a04, a02); \
  3005. MKL_DC_SUB_MUL_CONJ(a24, a24, a14, a12); \
  3006. MKL_DC_INV(a22, a22); \
  3007. MKL_DC_MUL(a23, a23, a22); \
  3008. MKL_DC_MUL(a24, a24, a22); \
  3009. dc_access(a, lda, 2, 3) = a23; \
  3010. dc_access(a, lda, 2, 4) = a24; \
  3011. MKL_DC_SUB_MUL_CONJ(a33, a33, a03, a03); \
  3012. MKL_DC_SUB_MUL_CONJ(a33, a33, a13, a13); \
  3013. MKL_DC_SUB_MUL_CONJ(a33, a33, a23, a23); \
  3014. if (MKL_DC_NON_POS(a33)) { \
  3015. dc_access(a, lda, 3, 3) = a33; \
  3016. *info = 4; \
  3017. break; \
  3018. } \
  3019. MKL_DC_SQRT(a33, a33); \
  3020. dc_access(a, lda, 3, 3) = a33; \
  3021. MKL_DC_SUB_MUL_CONJ(a34, a34, a04, a03); \
  3022. MKL_DC_SUB_MUL_CONJ(a34, a34, a14, a13); \
  3023. MKL_DC_SUB_MUL_CONJ(a34, a34, a24, a23); \
  3024. MKL_DC_INV(a33, a33); \
  3025. MKL_DC_MUL(a34, a34, a33); \
  3026. dc_access(a, lda, 3, 4) = a34; \
  3027. MKL_DC_SUB_MUL_CONJ(a44, a44, a04, a04); \
  3028. MKL_DC_SUB_MUL_CONJ(a44, a44, a14, a14); \
  3029. MKL_DC_SUB_MUL_CONJ(a44, a44, a24, a24); \
  3030. MKL_DC_SUB_MUL_CONJ(a44, a44, a34, a34); \
  3031. if (MKL_DC_NON_POS(a44)) { \
  3032. dc_access(a, lda, 4, 4) = a44; \
  3033. *info = 5; \
  3034. break; \
  3035. } \
  3036. MKL_DC_SQRT(a44, a44); \
  3037. dc_access(a, lda, 4, 4) = a44; \
  3038. } else { \
  3039. mkl_dc_type a10 = dc_access(a, lda, 1, 0); \
  3040. mkl_dc_type a20 = dc_access(a, lda, 2, 0); \
  3041. mkl_dc_type a30 = dc_access(a, lda, 3, 0); \
  3042. mkl_dc_type a40 = dc_access(a, lda, 4, 0); \
  3043. mkl_dc_type a21 = dc_access(a, lda, 2, 1); \
  3044. mkl_dc_type a43 = dc_access(a, lda, 4, 3); \
  3045. mkl_dc_type a31 = dc_access(a, lda, 3, 1); \
  3046. mkl_dc_type a41 = dc_access(a, lda, 4, 1); \
  3047. mkl_dc_type a32 = dc_access(a, lda, 3, 2); \
  3048. mkl_dc_type a42 = dc_access(a, lda, 4, 2); \
  3049. if (MKL_DC_NON_POS(a00)) { \
  3050. dc_access(a, lda, 0, 0) = a00; \
  3051. *info = 1; \
  3052. break; \
  3053. } \
  3054. MKL_DC_SQRT(a00, a00); \
  3055. dc_access(a, lda, 0, 0) = a00; \
  3056. MKL_DC_INV(a00, a00); \
  3057. MKL_DC_MUL(a10, a10, a00); \
  3058. MKL_DC_MUL(a20, a20, a00); \
  3059. MKL_DC_MUL(a30, a30, a00); \
  3060. MKL_DC_MUL(a40, a40, a00); \
  3061. dc_access(a, lda, 1, 0) = a10; \
  3062. dc_access(a, lda, 2, 0) = a20; \
  3063. dc_access(a, lda, 3, 0) = a30; \
  3064. dc_access(a, lda, 4, 0) = a40; \
  3065. MKL_DC_SUB_MUL_CONJ(a11, a11, a10, a10); \
  3066. if (MKL_DC_NON_POS(a11)) { \
  3067. dc_access(a, lda, 1, 1) = a11; \
  3068. *info = 2; \
  3069. break; \
  3070. } \
  3071. MKL_DC_SQRT(a11, a11); \
  3072. dc_access(a, lda, 1, 1) = a11; \
  3073. MKL_DC_SUB_MUL_CONJ(a21, a21, a20, a10); \
  3074. MKL_DC_SUB_MUL_CONJ(a31, a31, a30, a10); \
  3075. MKL_DC_SUB_MUL_CONJ(a41, a41, a40, a10); \
  3076. MKL_DC_INV(a11, a11); \
  3077. MKL_DC_MUL(a21, a21, a11); \
  3078. MKL_DC_MUL(a31, a31, a11); \
  3079. MKL_DC_MUL(a41, a41, a11); \
  3080. dc_access(a, lda, 2, 1) = a21; \
  3081. dc_access(a, lda, 3, 1) = a31; \
  3082. dc_access(a, lda, 4, 1) = a41; \
  3083. MKL_DC_SUB_MUL_CONJ(a22, a22, a20, a20); \
  3084. MKL_DC_SUB_MUL_CONJ(a22, a22, a21, a21); \
  3085. if (MKL_DC_NON_POS(a22)) { \
  3086. dc_access(a, lda, 2, 2) = a22; \
  3087. *info = 3; \
  3088. break; \
  3089. } \
  3090. MKL_DC_SQRT(a22, a22); \
  3091. dc_access(a, lda, 2, 2) = a22; \
  3092. MKL_DC_SUB_MUL_CONJ(a32, a32, a30, a20); \
  3093. MKL_DC_SUB_MUL_CONJ(a32, a32, a31, a21); \
  3094. MKL_DC_SUB_MUL_CONJ(a42, a42, a40, a20); \
  3095. MKL_DC_SUB_MUL_CONJ(a42, a42, a41, a21); \
  3096. MKL_DC_INV(a22, a22); \
  3097. MKL_DC_MUL(a32, a32, a22); \
  3098. MKL_DC_MUL(a42, a42, a22); \
  3099. dc_access(a, lda, 3, 2) = a32; \
  3100. dc_access(a, lda, 4, 2) = a42; \
  3101. MKL_DC_SUB_MUL_CONJ(a33, a33, a30, a30); \
  3102. MKL_DC_SUB_MUL_CONJ(a33, a33, a31, a31); \
  3103. MKL_DC_SUB_MUL_CONJ(a33, a33, a32, a32); \
  3104. if (MKL_DC_NON_POS(a33)) { \
  3105. dc_access(a, lda, 3, 3) = a33; \
  3106. *info = 4; \
  3107. break; \
  3108. } \
  3109. MKL_DC_SQRT(a33, a33); \
  3110. dc_access(a, lda, 3, 3) = a33; \
  3111. MKL_DC_SUB_MUL_CONJ(a43, a43, a40, a30); \
  3112. MKL_DC_SUB_MUL_CONJ(a43, a43, a41, a31); \
  3113. MKL_DC_SUB_MUL_CONJ(a43, a43, a42, a32); \
  3114. MKL_DC_INV(a33, a33); \
  3115. MKL_DC_MUL(a43, a43, a33); \
  3116. dc_access(a, lda, 4, 3) = a43; \
  3117. MKL_DC_SUB_MUL_CONJ(a44, a44, a40, a40); \
  3118. MKL_DC_SUB_MUL_CONJ(a44, a44, a41, a41); \
  3119. MKL_DC_SUB_MUL_CONJ(a44, a44, a42, a42); \
  3120. MKL_DC_SUB_MUL_CONJ(a44, a44, a43, a43); \
  3121. if (MKL_DC_NON_POS(a44)) { \
  3122. dc_access(a, lda, 4, 4) = a44; \
  3123. *info = 5; \
  3124. break; \
  3125. } \
  3126. MKL_DC_SQRT(a44, a44); \
  3127. dc_access(a, lda, 4, 4) = a44; \
  3128. } \
  3129. } while (0)
  3130. #define mkl_dc_potrf_uplo(uplo, n, a, lda, info, dc_access) \
  3131. do { \
  3132. if (n <= 5) { \
  3133. switch (n) { \
  3134. case 0: \
  3135. break; \
  3136. case 1: \
  3137. mkl_dc_potrf_1(uplo, n, a, lda, info, dc_access); \
  3138. break; \
  3139. case 2: \
  3140. mkl_dc_potrf_2(uplo, n, a, lda, info, dc_access); \
  3141. break; \
  3142. case 3: \
  3143. mkl_dc_potrf_3(uplo, n, a, lda, info, dc_access); \
  3144. break; \
  3145. case 4: \
  3146. mkl_dc_potrf_4(uplo, n, a, lda, info, dc_access); \
  3147. break; \
  3148. case 5: \
  3149. mkl_dc_potrf_5(uplo, n, a, lda, info, dc_access); \
  3150. break; \
  3151. } \
  3152. } else { \
  3153. mkl_dc_potrf_n(uplo, n, a, lda, info, dc_access); \
  3154. } \
  3155. } while (0)
  3156. #define mkl_dc_potrf_generic(uplo, n, a, lda, info, dc_access) \
  3157. do { \
  3158. int uuplo = MKL_DC_MisU(uplo); \
  3159. *info = 0; \
  3160. if (n != 0) { \
  3161. if (uuplo) { \
  3162. mkl_dc_potrf_uplo('U', n, a, lda, info, dc_access); \
  3163. } else { \
  3164. mkl_dc_potrf_uplo('L', n, a, lda, info, dc_access); \
  3165. } \
  3166. } \
  3167. } while (0)
  3168. static __inline void mkl_dc_potrf(char uplo, MKL_INT n, mkl_dc_type* a, MKL_INT lda, MKL_INT* info)
  3169. {
  3170. mkl_dc_potrf_generic(uplo, n, a, lda, info, MKL_DC_MN);
  3171. }
  3172. #ifndef MKL_DIRECT_CALL_LAPACKE_DISABLE
  3173. static __inline lapack_int mkl_dc_lapacke_potrf_convert(int matrix_layout, char uplo, lapack_int n, mkl_dc_type* a, lapack_int lda)
  3174. {
  3175. lapack_int info = 0;
  3176. if (MKL_DC_POTRF_CHECKSIZE(n)) {
  3177. if (matrix_layout == LAPACK_ROW_MAJOR) {
  3178. mkl_dc_potrf_generic(uplo, n, a, lda, &info, MKL_DC_MT);
  3179. } else {
  3180. mkl_dc_potrf_generic(uplo, n, a, lda, &info, MKL_DC_MN);
  3181. }
  3182. return info;
  3183. }
  3184. return MKL_DC_CONCAT3(LAPACKE_, MKL_DC_PREC_LETTER, potrf)(matrix_layout, uplo, n, a, lda);
  3185. }
  3186. #endif