#include <petsc.h>

int main(int argc, char *argv[])
{
   Mat A;
   Vec x, b;
   PetscErrorCode ierr;
   KSP ksp;
   PC pc;

   PetscInitialize(&argc, &argv, (char*)0, "void");
   
   ierr = MatCreate(PETSC_COMM_WORLD, &A); CHKERRXX(ierr);
   ierr = MatSetType(A,MATSEQAIJ); CHKERRXX(ierr);
   
   // Import matrix from file
   PetscViewer fileReader;
   ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD, "A", FILE_MODE_READ, &fileReader);
   CHKERRXX(ierr);
   ierr = MatLoad(A, fileReader); CHKERRXX(ierr);
   
   // KSP, PC
   ierr = KSPCreate(PETSC_COMM_WORLD, &ksp); CHKERRXX(ierr);
   ierr = KSPSetType(ksp, KSPGMRES); CHKERRXX(ierr);
   ierr = KSPSetOperators(ksp, A, A); CHKERRXX(ierr);
   ierr = KSPSetTolerances(ksp, 1e-6, 1e-9, PETSC_DEFAULT, 1000);
   CHKERRXX(ierr);
   ierr = KSPSetFromOptions(ksp); CHKERRXX(ierr);
   ierr = KSPGetPC(ksp, &pc); CHKERRXX(ierr);

   
   // Set up index sets for the Fieldsplit preconditioner
   PetscInt P_A_indices[324] = {224,226,228,230,232,234,236,238,240,242,244,246,248,250,252,254,256,258,260,262,264,266,268,270,272,274,276,278,280,282,284,286,288,290,292,294,328,330,332,334,336,338,340,342,344,346,348,350,352,354,356,358,360,362,364,366,368,370,372,374,376,378,380,382,384,386,388,390,392,394,396,398,432,434,436,438,440,442,444,446,448,450,452,454,456,458,460,462,464,466,468,470,472,474,476,478,480,482,484,486,488,490,492,494,496,498,500,502,536,538,540,542,544,546,548,550,552,554,556,558,560,562,564,566,568,570,572,574,576,578,580,582,584,586,588,590,592,594,596,598,600,602,604,606,640,642,644,646,648,650,652,654,656,658,660,662,664,666,668,670,672,674,676,678,680,682,684,686,688,690,692,694,696,698,700,702,704,706,708,710,744,746,748,750,752,754,756,758,760,762,764,766,768,770,772,774,776,778,780,782,784,786,788,790,792,794,796,798,800,802,804,806,808,810,812,814,848,850,852,854,856,858,860,862,864,866,868,870,872,874,876,878,880,882,884,886,888,890,892,894,896,898,900,902,904,906,908,910,912,914,916,918,952,954,956,958,960,962,964,966,968,970,972,974,976,978,980,982,984,986,988,990,992,994,996,998,1000,1002,1004,1006,1008,1010,1012,1014,1016,1018,1020,1022,1056,1058,1060,1062,1064,1066,1068,1070,1072,1074,1076,1078,1080,1082,1084,1086,1088,1090,1092,1094,1096,1098,1100,1102,1104,1106,1108,1110,1112,1114,1116,1118,1120,1122,1124,1126};
   PetscInt T_A_indices[324] = {225,227,229,231,233,235,237,239,241,243,245,247,249,251,253,255,257,259,261,263,265,267,269,271,273,275,277,279,281,283,285,287,289,291,293,295,329,331,333,335,337,339,341,343,345,347,349,351,353,355,357,359,361,363,365,367,369,371,373,375,377,379,381,383,385,387,389,391,393,395,397,399,433,435,437,439,441,443,445,447,449,451,453,455,457,459,461,463,465,467,469,471,473,475,477,479,481,483,485,487,489,491,493,495,497,499,501,503,537,539,541,543,545,547,549,551,553,555,557,559,561,563,565,567,569,571,573,575,577,579,581,583,585,587,589,591,593,595,597,599,601,603,605,607,641,643,645,647,649,651,653,655,657,659,661,663,665,667,669,671,673,675,677,679,681,683,685,687,689,691,693,695,697,699,701,703,705,707,709,711,745,747,749,751,753,755,757,759,761,763,765,767,769,771,773,775,777,779,781,783,785,787,789,791,793,795,797,799,801,803,805,807,809,811,813,815,849,851,853,855,857,859,861,863,865,867,869,871,873,875,877,879,881,883,885,887,889,891,893,895,897,899,901,903,905,907,909,911,913,915,917,919,953,955,957,959,961,963,965,967,969,971,973,975,977,979,981,983,985,987,989,991,993,995,997,999,1001,1003,1005,1007,1009,1011,1013,1015,1017,1019,1021,1023,1057,1059,1061,1063,1065,1067,1069,1071,1073,1075,1077,1079,1081,1083,1085,1087,1089,1091,1093,1095,1097,1099,1101,1103,1105,1107,1109,1111,1113,1115,1117,1119,1121,1123,1125,1127};
   PetscInt P_B_indices[352] = {0,2,4,6,8,10,12,14,16,18,20,22,24,26,28,30,32,34,36,38,40,42,44,46,48,50,52,54,56,58,60,62,64,66,68,70,72,74,76,78,80,82,84,86,88,90,92,94,96,98,100,102,104,106,108,110,112,114,116,118,120,122,124,126,128,130,132,134,136,138,140,142,144,146,148,150,152,154,156,158,160,162,164,166,168,170,172,174,176,178,180,182,184,186,188,190,192,194,196,198,200,202,204,206,208,210,212,214,216,218,220,222,296,298,300,302,304,306,308,310,312,314,316,318,320,322,324,326,400,402,404,406,408,410,412,414,416,418,420,422,424,426,428,430,504,506,508,510,512,514,516,518,520,522,524,526,528,530,532,534,608,610,612,614,616,618,620,622,624,626,628,630,632,634,636,638,712,714,716,718,720,722,724,726,728,730,732,734,736,738,740,742,816,818,820,822,824,826,828,830,832,834,836,838,840,842,844,846,920,922,924,926,928,930,932,934,936,938,940,942,944,946,948,950,1024,1026,1028,1030,1032,1034,1036,1038,1040,1042,1044,1046,1048,1050,1052,1054,1128,1130,1132,1134,1136,1138,1140,1142,1144,1146,1148,1150,1152,1154,1156,1158,1160,1162,1164,1166,1168,1170,1172,1174,1176,1178,1180,1182,1184,1186,1188,1190,1192,1194,1196,1198,1200,1202,1204,1206,1208,1210,1212,1214,1216,1218,1220,1222,1224,1226,1228,1230,1232,1234,1236,1238,1240,1242,1244,1246,1248,1250,1252,1254,1256,1258,1260,1262,1264,1266,1268,1270,1272,1274,1276,1278,1280,1282,1284,1286,1288,1290,1292,1294,1296,1298,1300,1302,1304,1306,1308,1310,1312,1314,1316,1318,1320,1322,1324,1326,1328,1330,1332,1334,1336,1338,1340,1342,1344,1346,1348,1350};
   PetscInt T_B_indices[352] = {1,3,5,7,9,11,13,15,17,19,21,23,25,27,29,31,33,35,37,39,41,43,45,47,49,51,53,55,57,59,61,63,65,67,69,71,73,75,77,79,81,83,85,87,89,91,93,95,97,99,101,103,105,107,109,111,113,115,117,119,121,123,125,127,129,131,133,135,137,139,141,143,145,147,149,151,153,155,157,159,161,163,165,167,169,171,173,175,177,179,181,183,185,187,189,191,193,195,197,199,201,203,205,207,209,211,213,215,217,219,221,223,297,299,301,303,305,307,309,311,313,315,317,319,321,323,325,327,401,403,405,407,409,411,413,415,417,419,421,423,425,427,429,431,505,507,509,511,513,515,517,519,521,523,525,527,529,531,533,535,609,611,613,615,617,619,621,623,625,627,629,631,633,635,637,639,713,715,717,719,721,723,725,727,729,731,733,735,737,739,741,743,817,819,821,823,825,827,829,831,833,835,837,839,841,843,845,847,921,923,925,927,929,931,933,935,937,939,941,943,945,947,949,951,1025,1027,1029,1031,1033,1035,1037,1039,1041,1043,1045,1047,1049,1051,1053,1055,1129,1131,1133,1135,1137,1139,1141,1143,1145,1147,1149,1151,1153,1155,1157,1159,1161,1163,1165,1167,1169,1171,1173,1175,1177,1179,1181,1183,1185,1187,1189,1191,1193,1195,1197,1199,1201,1203,1205,1207,1209,1211,1213,1215,1217,1219,1221,1223,1225,1227,1229,1231,1233,1235,1237,1239,1241,1243,1245,1247,1249,1251,1253,1255,1257,1259,1261,1263,1265,1267,1269,1271,1273,1275,1277,1279,1281,1283,1285,1287,1289,1291,1293,1295,1297,1299,1301,1303,1305,1307,1309,1311,1313,1315,1317,1319,1321,1323,1325,1327,1329,1331,1333,1335,1337,1339,1341,1343,1345,1347,1349,1351};

   IS P_A_IS, P_B_IS, T_A_IS, T_B_IS;
   ierr = ISCreateGeneral(PETSC_COMM_WORLD, 324,
                          P_A_indices, PETSC_COPY_VALUES, &P_A_IS);
   CHKERRXX(ierr);
   ierr = ISCreateGeneral(PETSC_COMM_WORLD, 324,
                          T_A_indices, PETSC_COPY_VALUES, &T_A_IS);
   CHKERRXX(ierr);
   ierr = ISCreateGeneral(PETSC_COMM_WORLD, 352,
                          P_B_indices, PETSC_COPY_VALUES, &P_B_IS);
   CHKERRXX(ierr);
   ierr = ISCreateGeneral(PETSC_COMM_WORLD, 352,
                          T_B_indices, PETSC_COPY_VALUES, &T_B_IS);
   CHKERRXX(ierr);
   
   ierr = PCFieldSplitSetIS(pc, "PA", P_A_IS); CHKERRXX(ierr);
   ierr = PCFieldSplitSetIS(pc, "TA", T_A_IS); CHKERRXX(ierr);
   ierr = PCFieldSplitSetIS(pc, "PB", P_B_IS); CHKERRXX(ierr);
   ierr = PCFieldSplitSetIS(pc, "TB", T_B_IS); CHKERRXX(ierr);
   ierr = KSPSetUp(ksp); CHKERRXX(ierr);

   VecCreate(PETSC_COMM_WORLD, &x);
   VecSetSizes(x, PETSC_DECIDE, 324*2 + 352*2);
   VecSetFromOptions(x);
   VecDuplicate(x, &b);
   VecSet(x, 1);
   MatMult(A,x,b);
   VecSet(x, 0);

   // Solve the system
   ierr = KSPSolve(ksp, b, x); CHKERRXX(ierr);

   PetscFinalize();
   return 0;
}